© The Authors 2025

1. Introduction

Weak lensing is the subtle distortion of light from distant galaxies through massive cosmic structures. This is a cornerstone for mapping the mass distribution in the Universe and for refining cosmological parameters. Upcoming cosmological surveys such as the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS), Hyper Suprime-Cam (HSC), Euclid, and Legacy Survey of Space and Time LSST (Heymans et al. 2012; Mandelbaum & Hyper Suprime-Cam HSC Collaboration 2017; Laureijs et al. 2011; Ivezić et al. 2019) aim to leverage weak lensing to address fundamental questions in cosmology (Lesgourgues & Pastor 2012; Li et al. 2019; Huterer 2010) and place tighter constraints on cosmological parameters (Troxel & Ishak 2015).

Power spectra are widely used tools for analysing weak gravitational lensing and extracting information about the matter distribution in the Universe. These second-order statistics have significant limitations, however: They fail to capture the non-Gaussian information generated by the non-linear evolution of structures on small scales, which leads to a substantial loss of insight into gravitational and astrophysical processes. An alternative is to use higher-order statistics, which can better capture the non-linear and non-Gaussian processes that govern cosmic structure formation. Examples include peak counts (Kruse & Schneider 1999; Liu et al. 2015a,b; Lin & Kilbinger 2015; Peel et al. 2017; Ajani et al. 2020), scattering transforms (Ajani et al. 2021; Cheng & Ménard 2021), Minkowski functionals (Kratochvil et al. 2012; Parroni et al. 2020), higher-order moments (Petri 2016; Peel et al. 2018; Gatti et al. 2020), wavelet ℓ₁-norms (Ajani et al. 2021; Sreekanth et al. 2024), one-point probability density functions (PDFs) (Bernardeau & Valageas 2000; Liu & Madhavacheril 2019; Barthelemy et al. 2021), and neural networks (Ribli et al. 2019; Fluri et al. 2018, 2019; Sharma et al. 2024; Lanzieri et al. 2023). These statistics enhance cosmological inference by capturing complex correlations and reducing degeneracies between cosmological parameters.

Accurately modelling these statistics, however, requires numerical simulations because analytical approaches often break down in non-linear regimes. Simulations are indeed the only practical means of capturing the complex and non-linear evolution of cosmological structures while incorporating the effects of baryons, biases, and astrophysical processes. They also enable the connection between theoretical predictions and real observations by accounting for systematic effects and uncertainties. Cosmological inference based on simulations is thus an essential step to maximise the information extracted from modern surveys and to test models of the Universe with precision. This has motivated different interesting new strategies in the past years for cosmological inferences, such as forward modelling (Alsing et al. 2016; Böhm et al. 2017; Remy et al. 2023; Lanzieri et al. 2023), field-level inference (Porqueres et al. 2023), or simulation-based inference (SBI; e.g. Gatti et al. 2024; Jeffrey et al. 2021; Papamakarios et al. 2017; Fagioli et al. 2020; Alsing et al. 2018; Tortorelli et al. 2020). In the forward-modelling approach, simulations are run to generate synthetic data based on a specific model and cosmological parameters. These synthetic datasets are then compared to real observational data. The comparison is often made using summary statistics (e.g. power spectra, correlation functions, or higher-order statistics), which are easier to handle computationally and can summarise the key features of the data.

The goal therefore is not to recover the full detailed matter density field, but rather to understand how theoretical models (with specific parameters) translate into measurable quantities that can be compared to observations. Field-level inference works directly with the raw data instead of relying on simplified summaries and can therefore recover the full matter density field. Forward-modelling and field-level inference both involve on-the-fly simulations. In SBI, simulations are usually pre-generated based on a range of cosmological models with different parameters (e.g. matter density or dark energy). These simulations provide synthetic datasets that reflect different possible configurations of the Universe according to the models being tested.

While on-the-fly simulations offer flexibility and adaptability for modelling complex systems in real time, their main limitations include high computational costs, potential loss of resolution and accuracy, challenges in exploring large parameter spaces efficiently, and difficulties in managing uncertainties. For SBI, its effectiveness is limited by the range and accuracy of the simulations it uses. If the simulations are not comprehensive or do not include all relevant physical processes, the inferences made could be biased or inaccurate. As the number of parameters increases, the complexity of generating and comparing simulations also grows significantly (i.e. the curse of dimensionality). Therefore, SBI limitations are mainly due to the high computational costs, the dependence on the quality of simulations, and challenges in exploring large parameter spaces.

This shows that these new ideas are very attractive, but still raise challenges, especially for the analysis of large stage-IV surveys. Emulators can serve as a middle ground between traditional methods based on second-order statistics and these recent simulation-intensive techniques. A huge variety of emulators exist in the literature (Feng et al. 2016; Izard et al. 2016; Tassev et al. 2013; White et al. 2014; Kitaura et al. 2014; Scoccimarro & Sheth 2002). Several approaches are grounded in analytical prescriptions that model the density field within redshift shells (Tessore et al. 2023) or employ inverse-Gaussianisation techniques (Yu et al. 2016). In recent years, several innovative methods have emerged that used machine-learning tools such as generative adversarial networks (GANs) (Boruah et al. 2024; Rodríguez et al. 2018) or normalising flows (Dai & Seljak 2022, 2023). They are trained using a set of pre-existing simulations. Nevertheless, the main limitation is that a large number of simulations are required for training, which constrains these methods to the parameter space in which they have been trained.

Another method for generating simulations of weak-lensing data relies on assuming a prior on the distribution of the convergence maps. The convergence map is indeed often described as a log-normal random field (Hilbert et al. 2011; Xavier et al. 2016). This can be very helpful for making more accurate inferences while reducing the reliance on full simulations, especially when the computational resources are limited.

We introduce a new way of generating weak-lensing convergence maps that can be obtained directly from the theory. This paper is organised as follows: We briefly review in Section 2 the current state of the weak-lensing emulators and overview the shifted log-normal models that are being used widely. In Section 3, we recapitulate the weak-lensing model, which is followed by the introduction of the algorithm in detail in Section 3.4. The full algorithmic procedure is outlined in Appendix A. In Section 4 we show the results from our algorithm and proceed to show the benchmarking in Section 4.3, before we finally conclude in Section 6. Appendix B provides complementary validation tests based on other higher-order statistics, namely one-point PDFs and peak counts. These were omitted from the main text for brevity. The appendix shows that these statistics also match remarkably well within the statistical errors of the numerical experiments.

2. Weak-lensing convergence map and the log-normal model

2.1. Convergence maps

The convergence κ(ϑ) represents the projection of the density field along the line of sight, weighted by a lensing kernel involving comoving distances (Mellier 1999),

$$\begin{array}{c}\hfill \kappa (\mathit{\vartheta })={\displaystyle {\int }_0^{{\chi }_s}}\mathrm{d}\chi \omega (\chi ,{\chi }_s)\delta (\chi ,\mathcal{D}\mathit{\vartheta }),\end{array}$$(1)

where χ is the comoving radial distance, χ_s is the radial distance to the source, and f_K(χ) is the comoving angular distance,

$$\begin{array}{c}\hfill f_K(\chi )=\{\begin{array}{cc}& \frac{sin(\sqrt{K}\chi )}{\sqrt{K}}\text{for}K>0,\hfill \\ \hfill & \chi \text{for}K=0,\hfill \\ \hfill & \frac{sinh(\sqrt{-K}\chi )}{\sqrt{-K}}\text{for}K<0,\hfill \end{array}\end{array}$$(2)

with K the spatial curvature. The lensing kernel ω(χ, χ_s) is defined as

$$\begin{array}{c}\hfill \omega (\chi ,{\chi }_s)=\frac{3{\mathrm{\Omega }}_m{H}_0^2}{2c^2}\frac{f_K(\chi )f_K({\chi }_s-\chi )}{f_K({\chi }_s)}(1+z(\chi )),\end{array}$$(3)

where c is the speed of light, Ω_m is the matter density parameter, and H₀ is the Hubble constant at z = 0.

The observed shear field γ(θ) relates to the convergence κ via mass inversion (Kaiser & Squires 1993; Starck et al. 2006; Martinet et al. 2015).

2.2. Angular Cl_s

Second-order statistics, such as the shear two-point correlation functions ξ_±(θ) and their Fourier-space counterpart, the angular power spectrum C_ℓ, play a crucial role in extracting cosmological information from weak-lensing surveys. These statistics describe the spatial distribution of the shear field and provide a window into the underlying matter distribution and its evolution.

Within the framework of the Limber approximation, the angular power spectrum of the convergence field for a specific tomographic bin is expressed as (Kilbinger 2015)

$$\begin{array}{c}\hfill C_{\kappa }(\ell )=\frac{9{\mathrm{\Omega }}_m^2{H}_0^4}{4c^4}{\displaystyle {\int }_0^{{\chi }_{\text{lim}}}}\text{d}\chi \frac{g^2(\chi )}{a^2(\chi )}{P}_{\delta }(k=\frac{\ell }{f_K(\chi )},\chi ),\end{array}$$(4)

where P_δ represents the matter power spectrum of the density contrast, χ denotes the comoving distance, and g(χ) is the lensing kernel that quantifies the lensing efficiency as a function of distance. The term a(χ) refers to the scale factor at a distance χ, while f_K(χ) is the comoving angular diameter distance, determined by the spatial curvature of the Universe. This formulation enables us to compute the convergence power spectrum, which provides insights into structure growth and the geometry of the Universe by linking weak-lensing observables to the underlying matter distribution.

2.3. Shifted log-normal model

In cosmology, random fields such as the matter density contrast or lensing convergence are also described statistically using probability density functions (PDFs), which define how likely different values of the field are. For some fields, the PDF is well approximated by a Gaussian distribution.

Many cosmological fields deviate significantly from this Gaussian behaviour, however. For instance, the matter density contrast and lensing convergence have hard lower limits (e.g. density cannot drop below zero relative to the average) and exhibit skewed shapes with heavy tails, indicating a higher likelihood of extreme values compared to a Gaussian distribution. In these cases, a shifted log-normal (Coles & Jones 1991; Xavier et al. 2016; Tessore et al. 2023) distribution offers a better model and has emerged as a practical alternative to full N-body simulations. This distribution can account for the asymmetry, hard limits, and heavy tails observed in the data, providing a more realistic representation of the statistical behaviour of these fields.

This approach is particularly effective in capturing the non-Gaussian distribution of matter while accounting for the negative convergence values observed in specific scenarios. By shifting the log-normal distribution, the model extends its range to include regions with negative convergence, thereby improving the realism and accuracy of density fluctuation simulations. This enhancement is particularly valuable for representing under-dense regions while preserving the non-linear clustering effects of large-scale structures.

Under moderate noise levels and typical accuracy requirements, the log-normal model provides reliable covariance matrices and bispectra (Hall & Taylor 2022). The skewness of the log-normal model alone is insufficient to fully capture the nonlinear features induced by gravitational evolution, however, leading to inaccuracies in higher-order statistics (Piras et al. 2023).

Given a set of variables K_gi following a multivariate Gaussian distribution with mean vector μ_i and covariance matrix ${\xi }_{\mathit{ij}}^g$, the shifted log-normal transformation is expressed as

$$\begin{array}{c}\hfill K_{log-n_i}=e^{K_{g_i}}-{\lambda }_i,\end{array}$$(5)

where λ_i depends on the cosmology and is computed as a function of the expected mean, variance, and skewness of the distribution. This dependence makes the model more flexible than a simple Gaussian approximation. Several approaches are used to fine-tune the shift parameter that is used in a log-normal model. The common approach is to focus on one scale of interest, and obtain a shift parameter such that the log-normal model has the correct skewness as that of the simulation, that is, it is calibrated against the PDF of the field at this particular scale, and that the shift parameter depends on the cosmology and on the scale. This means that fixing the shift parameter at one scale does not guarantee the exact skewness across all the other scales. Since convergence is not strictly a log-normal variable, different methods for determining λ_i yield different values, which might compromise the model accuracy at other scales. An added complexity arises when the matter density field is modelled as a log-normal field and integrated along the line of sight to derive the convergence map, which further complicates the estimation of λ_i. For further details, we refer to Xavier et al. (2016).

This log-normal transformation modifies the correlation function as follows (Xavier et al. 2016):

$$\begin{array}{c}\hfill {\xi }_{\mathit{ij}}^{\mathrm{logn}}={\lambda }_i(e^{{\xi }_{\mathit{ij}}^g}-1).\end{array}$$(6)

Hence, to generate a log-normal field with the desired correlation function, the Gaussian field to which the local transformation is applied must first be constructed with a corrected correlation function. As expected, Equation (6) shows that the correlation function (and hence, the variance) is independent of the shift parameter, which is often calibrated so as to reproduce the expected skewness of the log-normal field. The correlation function is related to the power spectrum via

$$\begin{array}{c}\hfill C_{\mathrm{logn}}^{\mathit{ij}}(\ell )=2\pi {\displaystyle {\int }_0^{\pi }}d\theta sin\theta P_{\ell }(cos\theta ){\xi }_{\mathit{ij}}^g(\theta ),\end{array}$$(7)

where P_ℓ represents the Legendre polynomial of order ℓ.

To illustrate the limitations inherent to the (shifted) log-normal model, we compared the convergence map generated using a log-normal model to a map derived from an N-body simulation. For this comparison, we used the publicly available Takahashi simulation suite (Takahashi et al. 2015), which provides full-sky convergence maps at a fixed cosmology¹. These maps are available at two grid resolutions of 4096 and 8192, and they include data spanning the full sky at intervals of 150 h⁻¹ comoving radial distance, covering redshifts from z = 0.05 to 5.3. The smoothing of the convergence maps is performed through a direct convolution of the map with a wavelet filter at specified angular scales. For this analysis, we used the convergence map at nside = 4096, which was subsequently downgraded to nside = 1024. A corresponding log-normal map with nside = 1024 and a maximum ℓ-mode of ℓ_max = 3 nside + 1 was then generated for comparison. Figure 1 compares the distribution of a log-normally generated convergence map with N-body simulations. The shift parameter was optimised to match the distribution of the N-body simulation at a smoothing scale of θ = 20 arcmin. The figure illustrates that this configuration fails to hold at other scales or for the full map, however.

Fig. 1. Comparison of the one-point PDF from N-body simulations and from maps generated by shifted log-normal modelling. Top panel: Distributions of the shifted log-normally generated maps (solid lines) and N-body simulations (dashed lines) for the unsmoothed case (in green) and at a smoothing angle of θ = 5 arcmin, θ = 12 arcmin, and at θ = 20 arcmin, as labelled. Bottom panel: Relative residuals of the log-normal PDF with respect to the N-body simulation.

This limitation underscores the necessity of developing a more flexible emulator capable of functioning consistently across all scales while effectively capturing cross-scale correlations. These advancements are crucial for robust multiscale analyses because they enable the extraction of non-linear information that is essential for cosmological inference.

3. Emulator based on the wavelet l1-norm theoretical prediction

Recent advancements in theoretical modelling have focused on higher-order statistics, such as the probability distribution function (PDF) (Barthelemy et al. 2021; Boyle et al. 2021) and the wavelet ℓ₁-norm (Sreekanth et al. 2024), which are based on large deviation theory (LDT). These models have demonstrated the capability to provide accurate predictions of the one-point statistics (PDF and wavelet ℓ₁-norm) of the density and convergence field up to mildly non-linear scales It was shown that for the wavelet ℓ₁-norm of the convergence maps, a prediction can be obtained at subpercent accuracy for a κ map smoothed up to θ = 15 arcmin for a source redshift z_s ≈ 2.

Our goal therefore was to build an emulator based on the theoretical prediction of the ℓ₁-norm. In the following section, we describe our algorithm in detail, which we call Generative mOdeLling of Convergence maps based ON predicteD one-point stAtistics (GOLCONDA), and we present our results.

3.1. Wavelet ℓ₁-norm

The wavelet transform is a powerful tool in the analysis of astronomical images because it can decompose data into components at different scales. This multi-scale decomposition is particularly effective for studying the complex hierarchical structures often found in astronomical data. We refer to Starck et al. (2015) for more details.

We employed a wavelet transform derived from using a tophat filter as the scaling function, as explained in (Sreekanth et al. 2024). An image c₀ can then be decomposed as the sum of all wavelet scales and a coarse resolution image c_J,

$$\begin{array}{c}\hfill c_0(x,y)=c_J(x,y)+{\displaystyle \sum _{j=1}^{J_{\text{max}}}}{w}_j(x,y),\end{array}$$(8)

where J_max denotes the maximum number of scales, and w_j represents wavelet images that capture the details of the original image at dyadic scales j. Each set of wavelet coefficients w_j was generated by convolving the input map with the respective wavelet kernel.

The wavelet ℓ₁-norm is then defined as

$$\begin{array}{c}\hfill {\ell }_1^{{\theta }_j}[i]={\displaystyle \sum _{k=1}^{\#coef(S_{{\theta }_j,i})}}|S_{{\theta }_j,i}[k]|,\end{array}$$(9)

where the set of coefficients at scale θ_j and amplitude bin i, S_θj, i = {w_θj / w_θj(x,y)∈[B_i,B_i + 1]}, depicts the wavelet coefficients w_θj having an amplitude within the bin [B_i, B_i + 1], and the pixel indices are given by (x, y). In other words, for each bin, the number of pixels k that fall in the bin waws collected, and the absolute values of these pixels were summed to obtain the ℓ₁-norm at this bin i.

This definition allowed us to capture data through the absolute values of all pixels within the map instead of limiting the analysis to identifying local minima or maxima. It offers the benefit of a multi-scale approach and was shown to be particularly robust compared to other approaches (Huber 1987; Giné et al. 2003).

3.2. Theoretical prediction of the wavelet ℓ₁-norm

Sreekanth et al. (2024) developed a theoretical framework to predict the wavelet ℓ₁-norm, denoted as ℓ₁^θ_j[i], for weak-lensing convergence maps. This framework leverages LDT, which is a statistical method for describing rare events, to derive the one-point probability density function (PDF) of the convergence field. In this framework, the wavelet ℓ₁-norm for a specific bin i is expressed as

$$\begin{array}{c}\hfill {\ell }_1^{{\theta }_j}[i]=|S_{{\theta }_j,i}|P(S_{{\theta }_j,i})\mathcal{N},\end{array}$$(10)

where |S_θj, i| is the bin value (absolute value of the wavelet coefficient), P(S_θj, i) is the normalised PDF of the wavelet coefficients at scale θ_j, and 𝒩 is a factor accounting for the normalisation of the predicted PDF.

The theoretical predictions for ℓ₁^θ_j[i] were validated against results from N-body simulations and demonstrated agreement at the percent level. This approach provides an efficient and accurate alternative to computationally expensive simulations for analysing non-Gaussian features in weak-lensing data. For further details, we refer to Sreekanth et al. (2024).

Given the simplicity of predicting the wavelet ℓ₁-norm for any specified cosmology, the objective was to use this prediction alongside the power spectrum as a constraint to generate a weak-lensing κ-map. The following section explores how these inputs can be effectively employed to produce a map that satisfies the desired statistical properties.

3.3. Optimisation problem

Our goal was to emulate a map with the statistical properties predicted by theoretical cosmological models for a given set of parameters. As the GLASS emulator, we wished to emulate a map κ that had the right theoretical power spectrum P_t, that is, 𝒫(κ) = P_t, where 𝒫 is the operator computing the radial power spectrum of κ. We also wished the same map to have the right distribution in the wavelet space as predicted by LDT. Noting, W_κ = {w₁,…,w_J − 1,c_J − 1} the wavelet transform of κ using J scales, that is, W = Φ^tκ, and $L_{1,t}=\{L_{1,t}^{{\theta }_1},\cdots ,L_{1,t}^{{\theta }_2},\dot{,}\cdots ,L_{1,t}^{{\theta }_J}\}$ the theoretical prediction in the different scales, each wavelet scale w_j should have the correct l₁-norm distribution ℓ₁^θ_j, i, that is, ℒ₁^θ_j(w_j) = L_1, t^θ_j, where ℒ₁^θ_j is the operator calculating the ℓ₁-norm at scale j following Eq. (9).

Specifically, we aimed to ensure that

• The map had the correct theoretical power spectrum P_t,

  $$\begin{array}{c}\hfill \mathcal{P}(\kappa )=P_t,\end{array}$$(11)

  where 𝒫 is the operator computing the radial power spectrum;

• The wavelet coefficients of the map had the correct ℓ₁-norm distribution at each scale,

  $$\begin{array}{c}\hfill \forall j\in \{1,\cdots ,J-1\},{\mathcal{L}}_1^{{\theta }_j}(w_j)=L_{1,t}^{{\theta }_j},\end{array}$$(12)

  where W_κ = {w₁, …, w_J − 1, c_J − 1}=Φ^tκ is the wavelet transform, and ℒ₁^θ_j is the operator computing the ℓ₁-norm at scale j.

We introduced two projection operators associated with these constraints:

• 𝒞_f: The projection onto the space of maps satisfying the power spectrum constraint (Eq. (11)).

• 𝒞_w: The projection onto the space of wavelet coefficients satisfying the ℓ₁-norm constraints (Eq. (12)).

The problem is then formulated as a constrained optimisation,

$$\begin{array}{c}\hfill \underset{\kappa ,W}{min}{\Vert \kappa -\mathrm{\Phi }W\Vert }^2\text{s.t.}{\mathcal{C}}_f(\kappa )=\kappa ,{\mathcal{C}}_w(W)=W.\end{array}$$(13)

This minimisation ensures that κ and W are consistent throughout the wavelet transform while satisfying the power spectrum and wavelet norm constraints. Because these constraints correspond to different domains (Fourier and wavelet space), a joint solution is not trivial and requires an iterative strategy. The algorithm called generalised forward-backward (GFB) (Liang et al. 2015) is well suited for problems involving multiple non-smooth convex constraints. It enabled us to decompose the optimisation into simpler sub-problems associated with each constraint and to alternate between them efficiently while ensuring convergence. In our case, the algorithm alternated between applying the projections 𝒞_f and 𝒞_w while enforcing consistency between κ and W.

An alternative minimisation approach works as follows:

• Assuming W is known, we minimise

  $$\begin{array}{c}\hfill \underset{\kappa }{min}\Vert \kappa -\mathrm{\Phi }W{\Vert }^{2}s.t.{\mathcal{C}}_f(\kappa )=\kappa .\end{array}$$(14)

• Assuming κ is known, we minimise

  $$\begin{array}{c}\hfill \underset{W}{min}\Vert \kappa -\mathrm{\Phi }W{\Vert }^{2}s.t.{\mathcal{C}}_w(W)=W.\end{array}$$(15)

The GFB algorithm applied to this setup yields the following iterative update scheme:

1. $z_1={\mathcal{C}}_p(\mathrm{\Phi }{W}^n)$

2. $z_2={\mathcal{C}}_w(W^n+{\mathrm{\Phi }}^t({\kappa }^n-\mathrm{\Phi }{W}^n))$

3. ${\kappa }^{n+1}=\frac{1}{2}(z_1+\mathrm{\Phi }{z}_2)$

4. $W^{n+1}=\frac{1}{2}({\mathrm{\Phi }}^t{z}_1+z_2)$.

Any additional summary statistics with theoretical prediction may be incorporated in a straightforward way by adding an additional constraint on κ.

The ${\mathcal{C}}_p$ operator consists of i) taking the Fourier transform $\widehat{\kappa }$ of κ, ii) computing its radial power spectrum P_κ, iii) multiplying the Fourier coefficients of $\widehat{\kappa }$ by $\sqrt{\frac{P_t[k]}{P_{\kappa }[k]}}$ (i.e. $z[u,v]=\sqrt{\frac{P_t[k]}{P_{\kappa }[k]}}\widehat{\kappa }[u,v]$, with $k=\sqrt{u^2+v^2}$), and iv) taking the inverse Fourier transform of z. Similarly, the ${\mathcal{C}}_w$ operator independently adjusts the ℓ₁-norm of all wavelet scales. The projection operators 𝒞_f and 𝒞_w are described in detail in Appendix A.

When Φ is orthogonal, we have Φ^tΦ = Id, and step 2 of the previous equation becomes

$$\begin{array}{cc}\hfill z_2& ={\mathcal{C}}_w(W^n+{\mathrm{\Phi }}^t{\kappa }^n-{\mathrm{\Phi }}^t\mathrm{\Phi }{W}^n)\hfill \\ \hfill & ={\mathcal{C}}_w({\mathrm{\Phi }}^t{\kappa }^n).\hfill \end{array}$$(16)

The method therefore becomes extremely simple by just iteratively applying the two steps

$$\begin{array}{cc}\hfill {\kappa }^{n+1}& =\frac{1}{2}({\mathcal{C}}_p(\mathrm{\Phi }{W}^n)+\mathrm{\Phi }{\mathcal{C}}_w({\mathrm{\Phi }}^t{\kappa }^n))\hfill \\ \hfill & =\frac{1}{2}({\mathcal{C}}_p({\kappa }^n)+\mathrm{\Phi }{\mathcal{C}}_w(W^n))\hfill \\ \hfill W^{n+1}& =\frac{1}{2}({\mathrm{\Phi }}^t{\mathcal{C}}_p(\mathrm{\Phi }{W}^n)+{\mathcal{C}}_w({\mathrm{\Phi }}^t{\kappa }^n))\hfill \\ \hfill & ={\mathrm{\Phi }}^t{\kappa }^{n+1},\hfill \end{array}$$(17)

and we fixed κ⁰ to a Gaussian or log-normal realisation, and W⁰ = Φ^tκ⁰. In the most general cases, the wavelet transform is a frame and is not necessarily orthogonal. We therefore have ΦΦ^t = Id (because we have a transform with an exact reconstruction), but we do not always have Φ^tΦ = Id. It remains a reasonable approximation, however, that allowed us to speed up the processing.

3.4. Algorithm

The algorithm iteratively adjusts an image to satisfy two main constraints: matching a target power spectrum through Fourier amplitude correction, and achieving desired ℓ₁-norm wavelet coefficients. This iterative process ensures that the final image accurately represents both the specified statistical properties and the target power spectrum.

To do this, the algorithm begins with an initial image, a realisation of a 2D Gaussian random field in our case. It proceeds by iteratively refining the image through a two-step process. First, an ℓ₁-norm-based wavelet coefficient adjustment (WCA) is applied to align the image wavelet coefficients with the desired properties. Second, the Fourier amplitudes are corrected to ensure that the image matches a specified target power spectrum. After each iteration, the refined image is obtained as a weighted average of the two corrected images. This iterative process continues until a convergence criterion, typically based on the relative root mean square (RMS) error between iterations, is satisfied or a predefined maximum number of iterations is reached. The algorithm involved here is given in Algorithm 1.

Fig. 2. Target convergence map and emulated convergence map obtained from our pipeline. The target map was obtained for a source redshift at zs = 2.134. The target and the generated maps both have a resolution of 1024 pixels on each side, with an angular size of 102 deg. The x- and y-axes give the number of pixels, and the pixel value is given by the color bar.

4. Results and validation of the emulator

4.1. Simulations

To test our pipeline, we used the suite called scinet light-cone simulations (SLICS) (Harnois-Déraps et al. 2018) of 924 fully independent N-body simulations that is publicly available². They were derived from a set of 1025 N-body simulations generated using the high-performance gravity solver CUBEP3M (Harnois-Déraps et al. 2013), which computed the nonlinear evolution of 1536³ particles starting from the initial conditions created at z = 120 using the Zeldovich approximation, in a box with a length of 505 h⁻¹Mpc on one side. The N-body code calculated the non-linear evolution of these collisionless particles up to z = 0. Using a multiple-plane technique, we generated convergence and shear maps over an area of 100 square degrees at 18 distinct source planes using the Born approximation.

4.2. Validation

To evaluate the performance of our emulator, we conducted a series of tests using a target image that was wavelet-decomposed into multiple scales employing a tophat filter as the scaling function. The set of equations governing the wavelet decomposition was previously defined in equation (8). As mentioned previously, our approach was motivated by recent theoretical advancements using the wavelet ℓ₁-norm, as predicted within the LDT framework (Sreekanth et al. 2024). These predictions provide input for generating emulated maps and are primarily valid in mildly non-linear scales, which consequently constrain the resolution of the maps. To demonstrate the effectiveness of this method, we validated our algorithm against the SLICS simulation suite.

We used a target image and decomposed it into a maximum of 2log(N) scales, where N is the number of pixels along one dimension of the map. This decomposition yielded 2log(N) sets of wavelet coefficients, along with a single coarse-scale image. The ℓ₁-norm values, computed per bin for both the wavelet coefficients and the coarse-scale image, were treated as the target wavelet distributions. Simultaneously, the C_ℓ values derived from the target image were designated as the target power spectra. Using these target distributions and power spectra as constraints, we applied the algorithm to generate emulated maps, which enabled us to assess its accuracy and robustness. For the validation, we generated a flat-sky map with a size of 100 deg² with 1024 pixels on each side.

Figure 2 illustrates the original target image compared to the emulated image generated through our pipeline. The target image was derived from the N-body simulation, while our software produced the reconstructed solution. At first glance, no significant visual differences are observed between the target and the reconstructed images up to the phases that are different in the two cases.

During the iterative process, the relative error in the wavelet ℓ₁-norms across different scales was assessed to evaluate the emulator performance. This error was calculated by comparing the ℓ₁-norm of the wavelet coefficients and the coarse-scale coefficient of the reconstructed image to those of the target image. Convergence was determined based on the stabilisation of these relative errors, indicating that the emulator had successfully approximated the target distributions.

The number of iterations in this work was fixed at n_iter = 150, as multiple tests indicated good convergence at this value. This conclusion was reached after verifying that the overall relative error, defined in Equation (18), in the ℓ₁-norm fell below 1%. This served as the convergence criterion to ensure that our iterative algorithm converged successfully. In our implementation, we set the tolerance Err_ℓ1 = 0.01, corresponding to a 1% relative error threshold. This value was chosen empirically, as we observed that both the ℓ₁-norm and the power spectrum residuals reliably stabilised below this level, indicating convergence of the key statistical properties. While tighter thresholds are possible, we found that Err_ℓ1 = 0.01 offered a practical trade-off between computational cost and statistical accuracy. The convergence behaviour is shown in Figure 3, where the relative error falls below 1% within approximately 100 iterations. Equation (18) computes the total normalised absolute error between the target and emulated wavelet ℓ₁-norms. For each scale j, it sums the absolute bin-wise differences over i, normalised by the maximum target value at that scale, and then sums over all scales.

$$\begin{array}{c}\hfill Er{r}_{{\ell }_1}={\displaystyle \sum _j}\left(\frac{{\displaystyle {\sum }_i}|{{\ell }_1}_{\text{target}}^{i,j}-{{\ell }_1}_{\text{emulated}}^{i,j}|}{\text{max}({{\ell }_1}_{\text{target}}^j)}\right).\end{array}$$(18)

Fig. 3. Evolution of the total relative error (in %) in the ℓ1-norm of the reconstructed emulated image compared to the simulation map during the iterative process. The total relative error proves a measure of how closely the reconstructed image approximates the target simulation map. The stabilisation of this error indicates the convergence of the emulator with the target distribution.

Figure 3 illustrates the evolution of the relative error over successive iterations. The plot highlights the error in the total ℓ₁-norm of the emulated image compared to the target image. This result demonstrates the convergence of the algorithm and its capability to replicate the target distributions at different scales with high accuracy. The plot shows that the algorithm converges quickly and reaches an accuracy at the percent level at around 100 iterations.

To conduct a more comprehensive comparison, we analysed the statistical properties of the reconstructed map, focusing on the power spectrum and the wavelet ℓ₁-norm. These comparisons are presented in Figures 4 and 5. Figure 4 presents the power spectra of the target κ map and the emulated map produced by our algorithm. The target power spectrum was derived from an N-body simulation, while the emulated power spectrum was generated using our pipeline, constrained by the input power spectrum and the wavelet ℓ₁-norm. The figure also shows the residuals of the power spectra, calculated as the relative difference between the target and emulated spectra across all angular multipoles. These residuals consistently remained at a sub-percent level, demonstrating that the algorithm can replicate the target power spectrum with high precision.

Fig. 4. Power spectra (Cℓ) of the target κ (solid orange line) map and the final emulated image (solid blue line) produced by our algorithm. The residuals (dashed green line) between the target and the emulated power spectra are also shown, indicating the level of agreement between the two maps.

Fig. 5. Wavelet ℓ1-norm at each scale and its residual between the decomposed target image (solid blue line) and the emulated image (dashed orange line) across different scales. The residuals, which are obtained as the difference between the target and the emulated wavelet ℓ1-norms, are also presented (dashed green line), reflecting the accuracy of the emulation process in reproducing the target distributions.

In Figure 5, we present the wavelet ℓ₁-norm of the decomposed image at various scales. The comparison highlights the wavelet ℓ₁-norm values computed for the target wavelet coefficients and those obtained from the emulated wavelet coefficients at each scale. The residuals, shown alongside, quantify the differences between the target and emulated wavelet ℓ₁-norms. Across all scales, the residuals are also at a sub-percent level, confirming the accuracy of the emulator in reproducing the target wavelet distributions.

These results emphasise the robustness of our algorithm, which successfully captures both the global and scale-dependent statistical properties of the target map. The sub-percent residuals observed in both the power spectrum and wavelet ℓ₁-norm validate the consistency of our approach and its potential for applications in cosmological analyses where precision is critical.

4.3. Benchmarking

In this section, we explore the time complexity, statistical properties, and higher-order statistics of the generated maps at multiple scales using different wavelet and smoothing filters (i.e. beyond the constraints used by the algorithm). This analysis provides insights into the strengths and limitations of the algorithm, highlighting its ability to emulate realistic data with minimum computational resources.

4.3.1. Performance

To generate the plots, our pipeline was able to generate each 500 × 500 map in less than 2 min. The Figure 6 below presents a plot of the processing time against map size. We used a single CPU and one core without incorporating any parallelisation or just-in-time (JIT) compilation, but it still performed comparably or better than existing methods.

Fig. 6. Time complexity of the generator. The x-axis shows the number of pixels on each side. On the y-axis, we plot the time taken to generate a map for the given size.

4.3.2. Higher cumulants

We further examined the statistical properties embedded in the maps, including skewness, variance, and kurtosis, at various scales over the entire map. To estimate the associated error bars, we divided each map into four patches and computed the standard deviation of the moments obtained across these patches. These results are displayed in Figure 7. The analysis was performed by dividing both the target and reconstructed maps into sub-patches, which allowed us to compute these statistics across the entire map and to estimate the associated uncertainties.

Fig. 7. Variance and skewness for the target image (in red) and the emulated image (in blue) obtained using a tophat scaling function to derive the wavelet. The emulated image was created by an algorithm that matched the wavelet ℓ1-norm of the target image from the SLICS simulation. Benchmarking was conducted by smoothing both images with different radii via a tophat filter, partitioning them into patches, and calculating skewness, variance, kurtosis, and their associated error bars by taking the standard deviation of the values obtained from four different patches of the maps. The x-axis represents the smoothing scale in arcminutes.

The left panel shows the skewness, which quantifies the asymmetry of the distribution of pixel values in the maps. The skewness values and associated errors of the target (red) and reconstructed (blue) maps agree closely, demonstrating the ability of the reconstruction process to preserve the asymmetry of the target map.

The middle panel presents the variance, which measures the dispersion of pixel values. As expected, the variance decreases with increasing smoothing scale due to the suppression of small-scale fluctuations during smoothing. The reconstructed map successfully reproduces the variance in the target map at all scales, with a good overlap between the red and blue points within the error bars.

The right panel shows the kurtosis, which captures the degree of tailedness in the distribution of the pixel values. The kurtosis generally decreases with increasing smoothing scale, consistent with the suppression of higher-order moments in smoothed data. The reconstructed map closely matches the kurtosis of the target map, with only minor deviations at some scales.

The results indicate that the reconstruction process preserves the higher-order statistical properties of the target map overall, as evidenced by the close agreement across skewness, variance, and kurtosis.

Appendix B presents some additional tests that omitted other higher-order statistics, namely PDFs and peak counts, showing that they incidentally match very well (within the statistical errors of the numerical experiment).

4.4. Extension to the generation of full-sky maps

Extending the algorithm to full-sky maps is a straightforward process because the method remains applicable in this setting. To do this, we used (Zonca et al. 2019; Górski et al. 2005). This extension has already been implemented, and the corresponding full-sky version is publicly available. Further details, including access to the implementation, can be found in Section ‘Data availability’.

5. Generalisation of the emulator to theoretical constraints

The validation presented in this work relies on a target wavelet ℓ₁-norm extracted from a simulation-generated weak- lensing convergence map. The method developed here is not restricted to this specific validation case, however. In principle, the emulator allows us to substitute the target values with any external input, including an arbitrary combination of a theoretical power spectrum C_ℓ and a corresponding set of ℓ₁-norm values across wavelet scales. This adaptability enables the generation of statistically controlled convergence maps, making it possible to impose user-defined constraints without an explicit dependence on simulation-based training data, which could come from purely theoretical grounds.

A notable advantage of this flexibility is the ability to incorporate stochasticity directly at the level of the input statistics so as to account for cosmic variance. In particular, any desired level of variance or non-Gaussian features arising from cosmic structure formation can be included by explicitly modulating the input target values. This approach allows us to control the sampling of cosmic variance, either by introducing uncertainty in the power spectrum C_ℓ or by perturbing the predicted ℓ₁-norms within physically motivated bounds. Consequently, the emulator can be extended to produce ensembles of statistically consistent weak-lensing maps that account for inherent fluctuations, thereby offering an alternative to expensive numerical simulations. One approach to achieve this might be to incorporate the cosmic variance at the level of the power spectrum only (see e.g. Moster et al. 2011) and obtain the corresponding wavelet ℓ₁-norm by only changing the power spectrum and keeping the scaled cumulants fixed Codis et al. (2016).

A fundamental limitation of using the theoretical prediction to generate maps is its dependence on scale. The derivation constrains the statistical properties of the wavelet ℓ₁-norm only up to mildly non-linear scales, as we discussed before. For scales smaller than the constrained range, however, the generated map may not accurately reproduce the statistical properties of the simulation.

Figure 8 shows an example of a map obtained by applying the constraints from the theory predictions at the first three scales (which includes the three wavelet coefficients and one coarse-scale coefficient). We also compared the wavelet ℓ₁-norm of the generated map with the theoretical predictions at the constrained scales, as shown in Figure 9. The figure clearly shows that the generated map successfully reproduces the theoretical ℓ₁-norm values with high accuracy. The residuals consistently remain at a sub-percent level, demonstrating the precision of the reconstruction and the effectiveness of the algorithm. The small oscillations in the residuals are attributed to the finite patch size used in the analysis.

Fig. 8. Emulated convergence map produced by our pipeline with a pixel size of 2 arcmin, generated by applying constraints on the first three scales for which the theory-based wavelet ℓ1-norm is obtained. The x- and y-axes represent the map range in degrees, and the pixel values are indicated by the colour bar.

Fig. 9. Comparison of the wavelet ℓ1-norm of the generated map (solid blue line) with the theoretical predictions (dashed orange line) across different scales ranging from 2 to 8 arcmin, shown from the top left to the bottom right panel. The final panel represents the coarse scale at 8 arcmin. The residuals (solid green line) between the target and emulated wavelet ℓ1-norms are also shown, indicating the accuracy of the emulation process in reproducing the expected distributions.

6. Conclusions

We developed and benchmarked a novel map-emulation algorithm that uses wavelet-based approaches, such as the tophat filter, to generate high-fidelity convergence maps with an enhanced accuracy and computational efficiency. Our method can emulate maps of 500 pixels per side in approximately one minute, even when it is restricted to a single CPU core and without parallelisation or just-in-time compilation. This performance underscores that our approach is more practical than traditional methods, which typically require substantial computational resources. The efficiency of the pipeline, coupled with its statistical rigour, provides a strong foundation for generating realistic cosmological maps in various applications.

A key focus of our work was the preservation and accurate reproduction of higher-order statistics, including variance, skewness, and kurtosis, at different scales and across entire maps. By conducting a thorough benchmarking analysis with a tophat filter, we demonstrated that our emulated maps closely approximate the statistical properties of target datasets, such as those derived from SLICS simulations. This capability to capture essential features of the underlying data makes our method particularly relevant for applications requiring precise modelling of small-scale structures and non-Gaussian features in cosmological fields.

Our results illustrate that the statistical alignment of the pipeline extends beyond basic moments, with a robust performance in reproducing higher-order metrics such as the cumulants across various smoothing scales apart from those used for constraining. This accuracy in representing complex structures suggests that our method provides a meaningful alternative to existing techniques, particularly log-normal map generation, which may not fully encapsulate higher-order dependencies and non-linear structures found in observed data.

Furthermore, the adaptability of our method offers promising opportunities for extensions and practical applications. The flexibility of the emulator enables the substitution of target values with any external input, including theoretical power spectra C_ℓ and corresponding sets of ℓ₁-norm values across wavelet scales. This adaptability allowed us to generate statistically controlled convergence maps with user-defined constraints without an explicit reliance on simulation-based training data. Additionally, the ability to incorporate stochasticity at the input level provides a means to model cosmic variance explicitly. This extends the ability of the emulator to produce ensembles of statistically consistent weak-lensing maps that account for inherent fluctuations.

While the current implementation of our algorithm is designed to reproduce specific one- and two-point statistics (namely, the wavelet ℓ₁-norm and the power spectrum) under idealised conditions, its modular structure allows a considerable flexibility. In particular, the framework can naturally be extended to incorporate observational systematics and physical effects. For instance, if theoretical predictions or simulation-based estimates of the wavelet ℓ₁-norm and power spectrum are available for convergence fields affected by intrinsic alignments or baryonic feedback, these can be directly used as constraints in the emulator. This enables the generation of convergence maps that encode the impact of these effects without requiring explicit modelling of the underlying processes. Moreover, although systematics such as shear calibration errors and photometric redshift uncertainties typically arise earlier in the weak-lensing analysis pipeline, the algorithm outputs can serve as inputs to downstream inference frameworks that incorporate these sources of uncertainty. Alternatively, the generated maps can be post-processed to simulate survey-specific features, including shape noise, masking, and selection biases. In summary, while this work has focused on validating the method in idealised scenarios, the design of the algorithm lends itself to future extensions that incorporate realistic systematics and survey conditions either through modified input constraints or post-processing.

Given these demonstrated capabilities, our wavelet-based emulation approach emerges as a robust and computationally efficient tool that might replace or complement log-normal map-generation techniques. By accurately preserving both higher-order statistics and fine-scale structures, our method offers a more complete representation of complex cosmological fields. Future developments may focus on optimising the algorithm through parallelisation and JIT techniques, expanding its usability and applicability across a wider range of datasets and scales. Furthermore, this method has potential applications in inpainting. Inpainting techniques aim to reconstruct missing data in an image or signal by estimating the absent values from the available information. Various approaches, such as sparse inpainting using the discrete cosine transform (Starck et al. 2005a,b), sparse recovery techniques (Pires et al. 2009; Lanusse et al. 2016), and Gaussian-constrained realisations (Zaroubi et al. 1995; Jeffrey et al. 2018), have been widely explored in the literature. The algorithm we presented might also prove useful for inpainting. While this exploration lies beyond the scope of this paper, assessing the effectiveness of this method compared to existing inpainting techniques remains an intriguing avenue for future research.

In conclusion, this work sets a foundation for efficient, accurate, and adaptable map-emulation strategies, offering significant potential for advancing cosmological data analysis, large-scale simulations, and other applications where statistical fidelity and computational efficiency are paramount.

Data availability

The code GOLCONDA (Generative mOd-eLing of Convergence maps based ON predicteD one-point stAtistics) is made available in the spirit of open research. It enabled us to reproduce the plots and can be accessed at GitHub.

Acknowledgments

This work was supported by the TITAN ERA Chair project (contract no. 101086741) within the Horizon Europe Framework Program of the European Commission, and the Agence Nationale de la Recherche (ANR-22-CE31-0014-01 TOSCA). We would like to thank Joachim Harnois-Déraps for making public the SLICS mock data, which can be found at http://slics.roe.ac.uk/. We would also like to thank Takahashi and collaborators for making public the TAKAHASHI simulation suite, which can be found at http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allsky_raytracing/. The authors thank Andreas Tersenov, Cora Uhlemann and Alexandre Barthelemy for useful discussions.

References

Appendix A: Algorithm in detail

A.1. ℓ₁-Norm matching wavelet coefficient adjustment algorithm

We are using the ℓ₁-Norm Matching Wavelet Coefficient Adjustment Algorithm (ℓ₁-WCA), which is designed to adjust wavelet coefficients at various scales to achieve specified ℓ₁ norms within predefined histogram bins. This algorithm iteratively modifies the wavelet coefficients in each bin, separately adjusting positive and negative values to meet target ℓ₁ norms while ensuring the values remain within set boundaries. We believe it could be useful in applications where specific statistical properties of wavelet-transformed data must be preserved, providing a flexible method to enforce norm constraints without disrupting the underlying data distribution.

Given a target L_1, t, ℓ₁ of the wavelet coefficient can be adjusted to match the target as follows:

First, the total adjustment Δ_total is obtained as the difference between the target and the current wavelet ℓ₁-norm.

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\text{total}}& =L_{1,t}^{{\theta }_j}-{\displaystyle \sum _{i=1}^{N_{\text{bin}}}}|w_{j,i}|,\hfill \end{array}$$(A.1)

where L_1, t^θ_j represents the target ℓ₁-norm at scale j, and the summation computes the current ℓ₁-norm of the wavelet coefficients. This is used to get the adjustment per positive coefficient as:

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_p& =\frac{{\mathrm{\Delta }}_{\text{total}}}{N_p},\hfill \end{array}$$(A.2)

where N_p is the number of positive coefficients. The adjusted positive coefficients are then given by:

$$\begin{array}{cc}\hfill w_{j,i}& =max(w_{j,i}+{\mathrm{\Delta }}_p,0),\forall i\text{such that}{w}_{j,i}\ge 0.\hfill \end{array}$$(A.3)

After this update, the new ℓ₁ is computed as:

$$\begin{array}{cc}\hfill L_{1,\text{updated}}& ={\displaystyle \sum _{i=1}^{N_{\text{bin}}}}|w_{j,i}|.\hfill \end{array}$$(A.4)

The remaining adjustment to the target norm is obtained as:

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\text{residual}}& =L_{1,t}^{{\theta }_j}-L_{1,\text{updated}}.\hfill \end{array}$$(A.5)

The residual is distributed amongst the negative coefficients and updated as:

$$\begin{array}{cc}\hfill w_{j,i}& =min(w_{j,i}-{\mathrm{\Delta }}_n,0),\forall i\text{such that}{w}_{j,i}<0.\hfill \end{array}$$(A.6)

This now gives a wavelet coefficient with the constrained ℓ₁-norm matching the target.

The corresponding algorithm is given in Algorithm A.1.

A.2. Algorithm for getting the correct power spectra

To generate an image with the desired target power spectrum (Cls), we perform Fourier amplitude rescaling. This process adjusts the Fourier amplitudes of the image to match the target power spectrum while preserving the original phase information to maintain spatial structure, as shown in Algorithm 3.

Appendix B: Evaluation of higher-order statistics

We conducted a detailed comparison of the emulated convergence maps with the target maps, focusing on the one-point Probability Distribution Function (PDF) and peak counts. In addition to the main analysis presented in the main text, we have also conducted a detailed comparison of the emulated convergence maps with the target maps using other high-order statistics. We have especially focused on the one-point Probability Distribution Function (PDF) and peak counts, which are detailed below.

B.1. Probability distribution function (PDF)

Figure B.1 shows the one-point PDF of the convergence field at two different smoothing scales (2 arcmin and 12 arcmin), obtained using a circular tophat filter. These filters are distinct from the wavelet-like difference-of-top hats used in defining the ℓ₁-norm constraint. The top panels show the PDFs for both the target and emulated maps, while the bottom panels display the residuals. The results indicate agreement at the ∼10% level across most of the distribution, with slight discrepancies in the extreme tails. Let us note that for intermediate to large smoothing scales (e.g. 12 arcmin here), the difference is entirely within the statistical error bars of our measurements, hence no significant departure is observed in this regime. For smaller smoothing scales, the deviation is more significant, but this is in a regime where non-linearities start to be quite important (hence less interesting for cosmological analysis if systematics can not be controlled), and resolution effects might also play a significant role.

Fig. B.1. Top panels: One-point PDF of the target and emulated maps at two smoothing scales (2 and 12 arcmin), computed using a circular tophat filter. Bottom panels: Residuals between emulated and target PDFs.

B.2. Peak count statistics

In Figure B.2, we show the peak count distribution at various wavelet scales. Peaks are defined as local maxima in the map, and the histograms are shown separately for the target and emulated maps. The agreement is very good for all scales, with some fluctuation at the largest scales (4 and 5), mostly driven by a lack of statistics.

Fig. B.2. Peak count distribution at different wavelet scales as labelled for the target (dotted red line) and emulated maps (blue solid line). Emulated maps reproduce the peak counts of the target with high fidelity at small to intermediate scales.

B.3. Implications for higher-order moments

Although the algorithm is optimised only to match the power spectrum and wavelet ℓ₁-norm, the consistency seen in the PDFs and peak counts, especially in the tails, indicates that higher-order moments are approximately preserved. These statistics, especially the peak counts, are sensitive to non-Gaussian features and implicitly reflect the presence of higher-order correlations beyond the fourth moment. However, we emphasise that our method does not guarantee fidelity in any specific higher-order moment. The observed agreement suggests the approach generalises well and captures broader statistical features, but also leaves room for future improvements that target additional moments explicitly.

All Figures

Fig. 1. Comparison of the one-point PDF from N-body simulations and from maps generated by shifted log-normal modelling. Top panel: Distributions of the shifted log-normally generated maps (solid lines) and N-body simulations (dashed lines) for the unsmoothed case (in green) and at a smoothing angle of θ = 5 arcmin, θ = 12 arcmin, and at θ = 20 arcmin, as labelled. Bottom panel: Relative residuals of the log-normal PDF with respect to the N-body simulation.

In the text

	Fig. 2. Target convergence map and emulated convergence map obtained from our pipeline. The target map was obtained for a source redshift at zs = 2.134. The target and the generated maps both have a resolution of 1024 pixels on each side, with an angular size of 102 deg. The x- and y-axes give the number of pixels, and the pixel value is given by the color bar.
In the text

	Fig. 3. Evolution of the total relative error (in %) in the ℓ1-norm of the reconstructed emulated image compared to the simulation map during the iterative process. The total relative error proves a measure of how closely the reconstructed image approximates the target simulation map. The stabilisation of this error indicates the convergence of the emulator with the target distribution.
In the text

	Fig. 4. Power spectra (Cℓ) of the target κ (solid orange line) map and the final emulated image (solid blue line) produced by our algorithm. The residuals (dashed green line) between the target and the emulated power spectra are also shown, indicating the level of agreement between the two maps.
In the text

Fig. 5. Wavelet ℓ1-norm at each scale and its residual between the decomposed target image (solid blue line) and the emulated image (dashed orange line) across different scales. The residuals, which are obtained as the difference between the target and the emulated wavelet ℓ1-norms, are also presented (dashed green line), reflecting the accuracy of the emulation process in reproducing the target distributions.

In the text

	Fig. 6. Time complexity of the generator. The x-axis shows the number of pixels on each side. On the y-axis, we plot the time taken to generate a map for the given size.
In the text

Fig. 7. Variance and skewness for the target image (in red) and the emulated image (in blue) obtained using a tophat scaling function to derive the wavelet. The emulated image was created by an algorithm that matched the wavelet ℓ1-norm of the target image from the SLICS simulation. Benchmarking was conducted by smoothing both images with different radii via a tophat filter, partitioning them into patches, and calculating skewness, variance, kurtosis, and their associated error bars by taking the standard deviation of the values obtained from four different patches of the maps. The x-axis represents the smoothing scale in arcminutes.

In the text

	Fig. 8. Emulated convergence map produced by our pipeline with a pixel size of 2 arcmin, generated by applying constraints on the first three scales for which the theory-based wavelet ℓ1-norm is obtained. The x- and y-axes represent the map range in degrees, and the pixel values are indicated by the colour bar.
In the text

Fig. 9. Comparison of the wavelet ℓ1-norm of the generated map (solid blue line) with the theoretical predictions (dashed orange line) across different scales ranging from 2 to 8 arcmin, shown from the top left to the bottom right panel. The final panel represents the coarse scale at 8 arcmin. The residuals (solid green line) between the target and emulated wavelet ℓ1-norms are also shown, indicating the accuracy of the emulation process in reproducing the expected distributions.

In the text

	Fig. B.1. Top panels: One-point PDF of the target and emulated maps at two smoothing scales (2 and 12 arcmin), computed using a circular tophat filter. Bottom panels: Residuals between emulated and target PDFs.
In the text

	Fig. B.2. Peak count distribution at different wavelet scales as labelled for the target (dotted red line) and emulated maps (blue solid line). Emulated maps reproduce the peak counts of the target with high fidelity at small to intermediate scales.
In the text