© The Authors 2026

1. Introduction

The Epoch of Reionisation (EoR) is a critical transition in the history of the Universe, marking the emergence of the first stars and the eventual ionisation of the once-neutral intergalactic medium (IGM). Indirect observations of the EoR from the cosmic microwave background (CMB) (Planck Collaboration XLVII 2016; Gorce et al. 2018) and the Gunn-Peterson trough in quasar spectra (Fan et al. 2006; Becker et al. 2015; Bosman et al. 2018; Qin et al. 2021), constrain the timing of the EoR between redshifts 14 ≳ z ≳ 5.3.

Ground–based CMB measurements provide a complementary EoR probe: South Pole Telescope (SPT) analyses constrain the duration to $\mathrm{\Delta }{z}_{\mathrm{re},50}=1.1_{-0.7}^{+1.6}$ (< 4.1 at 95% C.L.), and a cosmology–dependent reanalysis finds Δz_re, 50 < 4.3 with $z_{\mathrm{re}}=7.9_{-1.3}^{+1.1}$ (Reichardt et al. 2021; Gorce et al. 2022); although the Atacama Cosmology Telescope (ACT) has not yet detected the trispectrum, it sets an upper limit and roadmap for future detections (MacCrann et al. 2024), reinforcing the high-ℓ trispectrum as a viable, complementary EoR probe (Kumar et al. 2025).

While these indirect observations offer loose timings of the EoR, they do not constrain the sources responsible for the EoR, nor the state of the IGM, for which a direct probe is needed. The redshifted 21 cm signal, resulting from the hyperfine transition in neutral hydrogen during the EoR, presents a promising direct observational probe (observed by radio interferometers between ∼100 MHz and 200 MHz) for a comprehensive study of this epoch.

Observations of the 21 cm signal (EoR signal, henceforth) by radio interferometers must contend with foregrounds such as synchrotron emission from diffuse clouds in the Galaxy, noise arising from the antennas, and systematics such as mode-mixing arising from the chromatic nature of radio interferometers. All of these contaminants are orders of magnitude stronger in total variance than the signal itself (Shaver et al. 1999; Jelić et al. 2008; Hothi et al. 2021).

The cylindrically averaged power spectrum of a radio interferometer observation (Fig. 1) has three distinct regions in the k_∥ − k_⊥ plane¹ (Datta et al. 2010; Vedantham et al. 2012). The first region, located at low k_∥, is dominated by astrophysical foregrounds. These foregrounds include Galactic emission, arising from free-free bremsstrahlung and diffuse synchrotron emission, as well as extragalactic synchrotron emission arising from sources such as radio haloes and clusters (Jelić et al. 2008; Hothi et al. 2021). The second region is the foreground wedge, which arises from the chromatic nature of the instrument: the interferometer’s point spread function (PSF) is frequency dependent, leading to mode mixing between k_⊥ and k_∥ and causing the foreground emission to leak to higher k_∥ (Liu et al. 2014). Bright sources entering the side lobes of the interferometer’s PSF also contribute to the emission within the wedge, all of which are bound by the horizon². The third region is the EoR window, where foregrounds and instrumental effects are subdominant (there is little to no contribution on these scales) to the noise and EoR signal (see Munshi et al. 2025 for a more comprehensive overview of the wedge).

The EoR signal is expected to be highly non-Gaussian, and this non-Gaussian information is important to constrain both the astrophysics of this epoch and the progression of reionisation (Shimabukuro et al. 2017; Watkinson et al. 2022; Tiwari et al. 2022). To fully characterise the signal requires statistics beyond the power spectrum.

Scattering transforms (STs) have emerged as a novel set of summary statistics constructed by successive wavelet convolutions and the application of non-linear operators such as modulus (Mallat 2011; Bruna & Mallat 2012). They have emerged as powerful low-variance statistics to characterise non-Gaussian processes and have already been successfully applied to various topics, including the study of the interstellar medium (Allys et al. 2019; Regaldo-Saint Blancard et al. 2020; Lei & Clark 2023); weak lensing (Cheng et al. 2020); the large-scale structure of the Universe for 2D (Allys et al. 2020) and 3D data (Eickenberg et al. 2022; Valogiannis & Dvorkin 2022; Régaldo-Saint Blancard et al. 2024); and more recently, to the EoR (Greig et al. 2022; Zhao et al. 2024; Hothi et al. 2024; Prelogović & Mesinger 2024; Shimabukuro et al. 2025).

An advantage of ST statistics is that they can be used to build generative models of physical processes within the framework of maximum entropy models. Such models, parametrised by the ST statistics themselves, efficiently approximate complex non-Gaussian physical fields, even when constructed from a single sample (Allys et al. 2020; Cheng et al. 2023; Mousset et al. 2024; Campeti et al. 2025). Scattering transforms (STs) have also been used in the context of EoR to build diffusion models, using U-Nets, and outperform generative adversarial networks when applied to 2D EoR fields (Zhao et al. 2023).

The ability to generate realistic realisations of various processes under ST constraints has led to the development of new component separation algorithms. These ST-based component separations have already been successfully applied directly to observational data, exploiting the non-Gaussian properties of the components within the data, without physically driven priors on the component of interest (Regaldo-Saint Blancard et al. 2021; Delouis et al. 2022; Auclair et al. 2024). However, these previous algorithms have only been applied to 2D data. For EoR studies, this formalism must be extended to 3D EoR lightcones. This in turn requires reliance on ST statistics that can adequately describe such data, the development of which is the purpose of this paper. This work therefore extends ST generative models to 3D EoR data and quantitatively validates their quality.

To construct such 3D ST statistics, the first step is to introduce a relevant set of wavelets that characterise the different regions of the k_∥ − k_⊥ plane (see Fig. 1). In previous 3D ST applications, wavelet sets typically probe extended regions of the k_∥ − k_⊥ plane, thus mixing regions where the foregrounds are dominant and subdominant (see, for example, Eickenberg et al. 2022). This is problematic because mixing regions where foregrounds are dominant with subdominant regions is detrimental to the efficient extraction of information from the EoR signal (Blamart & Liu 2025). To avoid this drawback, we constructed a new set of 3D wavelets from the tensor product of 1D wavelets along the k_∥ axis and 2D isotropic wavelets along the sky axes (k_x, k_y). We then used this new wavelet set to construct 3D ST statistics, which we in turn used to construct a generative model of EoR lightcones, the quality of which was assessed quantitatively. The ST model is built from a single simulated EoR lightcone; we compared its statistical properties with those of the simulation using several independent and complementary statistics.

Fig. 1. Schematic of the cylindrically-averaged power spectrum of an observed EoR lightcone, with the three distinct regions highlighted.

The paper is structured as follows. Section 2 outlines the simulations used. Section 3 presents the 3D scattering transforms. Section 4 describes the generative model built from scattering transforms. Section 5 compares the synthesised lightcones with the target EoR lightcones.

2. LoReLi II simulation

The simulation used in this work is from the LoReLi II (LOw REsolution LIcorice) database (Meriot & Semelin 2024; Meriot et al. 2025). This database is a collection of Licorice simulations designed to study the EoR signal (Semelin et al. 2007, 2017).

The LoReLi II database comprises 9828 simulations of the EoR signal brightness temperature and explores a four-parameter space: the gas-to-star conversion timescale, the minimum halo mass for star formation, the escape fraction of UV radiation, and the X-ray production efficiency, while keeping the cosmological parameters constant at {H₀: 67.8 km/s/Mpc, Ω₀: 0.308, Ω_b: 0.0484, Ω_Λ: 0.692, and σ₈: 0.815}. This setup allows for a detailed exploration of the astrophysical processes during reionisation without the confounding effects of varying cosmology.

In this paper, the ST generative model is built from a single LoReLi simulation of a lightcone. The astrophysical parameters of the lightcone, which we refer to as the target lightcone, are the log10 of the minimum halo mass, M_min = 8.0; the X-ray luminosity fraction, f_x = 0.1; the hard X ray fraction, $r_{\frac{H}{S}}$ = 0.2; the gas conversion timescale for star formation in units of 10 Myr, τ = 8191.9; and the ionising escape fraction, f_esc, post = 0.275.

This lightcone spans a redshift range of 23 ≤ z ≤ 5. The xy-plane (the 2D sky) contains 256² pixels, and the z-domain contains 3072 redshift channels between z = 23 and z = 5. To reduce the computational cost and limit the z evolution over the lightcone, we downgraded it from 256² pixels to 128² pixels in the xy plane, and restricted the z values to 9.5 ≤ z ≤ 8.5 (frequency bandwidth δν ∼ 14 MHz), corresponding to 170 pixels in the z domain, for a final data cube of (128, 128, 170) pixels. If the ST statistics were estimated on the initial lightcone, their statistical properties would evolve with z, whereas we approximated only limited z-evolution in the final constructed lightcone. However, we note that a lightcone with a larger redshift range could be built by stacking models constructed with a more limited redshift range (Blamart & Liu 2025).

3. Scattering transforms for 3D lightcones

3.1. Wavelets for 3D lightcones

In this paper, we introduce a wavelet set that characterises the k_∥ − k_⊥ plane by performing a tensor product between wavelet sets defined on different axes. The first is a 1D wavelet set defined in k_z, denoted by ∥, and the second is a 2D isotropic wavelet set defined in the (k_x, k_y) plane, denoted by ⊥. To simplify notation, we define the modulus wavenumbers in the two _∥ and _⊥ directions by k_∥ = |k_z| and $k_{\perp }=\sqrt{k_x^2+k_y^2}$.

The L1-normalised mother wavelet of the 1D ∥ wavelet set is defined as

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_0^{\Vert }(k_z)={\psi }_0^{\Vert }(k_{\Vert })=\frac{1}{\sqrt{2\pi {\sigma }_0^{\Vert 2}}}exp-\frac{{(k_{\Vert }-{\mu }_0^{\Vert })}^2}{2{\sigma }_0^{\Vert 2}},\end{array}$$(1)

where μ₀^∥ is the central wavenumber of the mother wavelet in Fourier space, and σ₀^∥ its width. Here, we denote the wavelets in the 3D (k_x, k_y, k_z) space with Ψ and the wavelets in the 2D (k_∥, k_⊥) cylindrical plane with ψ. The L1-normalised mother wavelet of the 2D ⊥ isotropic wavelet set is defined as

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_0^{\perp }(k_x,k_y)={\psi }_0^{\perp }(k_{\perp })=\frac{1}{2\pi {\sigma }_0^{\perp 2}}exp-\frac{{(k_{\perp }-{\mu }_0^{\perp })}^2}{2{\sigma }_0^{\perp 2}},\end{array}$$(2)

where μ₀^⊥ is again the central wavenumber of the mother wavelet in Fourier space, and σ₀^⊥ its width.

We then constructed the two wavelet sets by dilating the respective mother wavelets,

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{j_{\Vert }}^{\Vert }(k_z)={\psi }_{j_{\Vert }}^{\Vert }(k_{\Vert })& =2^{-\frac{1}{2}\frac{j_{\Vert }}{Q_{\Vert }}}{\psi }_0^{\Vert }\left(2^{-\frac{j_{\Vert }}{Q{j}_{\Vert }}}{k}_{\Vert }\right),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{j_{\perp }}^{\perp }(k_x,k_y)={\psi }_{j_{\perp }}^{\perp }(k_{\perp })& =2^{-\frac{1}{2}\frac{2j_{\perp }}{Q{j}_{\perp }}}{\psi }_0^{\perp }\left(2^{-\frac{j_{\perp }}{Q_{\perp }}}{k}_{\perp }\right).\hfill \end{array}$$(4)

Here, j_∥ and j_⊥ are the integer scales of the dilation, and Q_∥ and Q_⊥ are the global quality factors. From these equations, ψ_j∥^∥(k_z) and ψ_j⊥^⊥(k_⊥) can be written in the same form as the mother wavelets defined in Eqs. (1) and (2) by making the following replacements:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\sigma }_0^{\alpha }\to {\sigma }_{j_{\alpha }}^{\alpha }=2^{\frac{j_{\alpha }}{Q_{\alpha }}}{\sigma }_0^{\alpha },\\ \hfill {\mu }_0^{\alpha }\to {\mu }_{j_{\alpha }}^{\alpha }=2^{\frac{j_{\alpha }}{Q_{\alpha }}}{\mu }_0^{\alpha },\end{array}\end{array}$$(5)

for α ∈ { ∥ , ⊥ }. The j_α integer scales of the ψ_jα wavelets range from 0 to J_α − 1, where J_α is the number of scales. Once J_α is chosen, the value of Q_α is fixed so that the wavelet set samples the full extent of the Fourier space. Specifically, when the mother wavelet, centred at μ₀^α, is dilated up to the largest scale, the resulting ψ_{Jα − 1} wavelet is centred at the maximum Fourier scale k_α, max, i.e. the Nyquist scale, yielding

$$\begin{array}{c}\hfill log\frac{k_{\alpha ,\max }}{{\mu }_0^{\alpha }}=\frac{J_{\alpha }-1}{Q_{\alpha }}log2.\end{array}$$(6)

To construct the 3D wavelet, we consider the tensor product between the two wavelet sets we previously defined:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{j_{\Vert },j_{\perp }}(k_x,k_y,k_z)={\mathrm{\Psi }}_{j_{\Vert }}^{\Vert }(k_z)& \otimes {\mathrm{\Psi }}_{j_{\perp }}^{\perp }(k_x,k_y)\hfill \\ \hfill & ={\psi }_{j_{\Vert },j_{\perp }}(k_{\Vert },k_{\perp })={\psi }_{j_{\Vert }}^{\Vert }(k_{\Vert })·{\psi }_{j_{\perp }}^{\perp }(k_{\perp }),\hfill \end{array}$$(7)

where ⊗ is the tensor product. In the following, we use the shorthand notation λ = {j_∥, j_⊥}, which yields

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\lambda }(k_x,k_y,k_z)={\psi }_{\lambda }(k_{\Vert },k_{\perp }).\end{array}$$(8)

These wavelets are real in Fourier space and complex in real space.

In this paper, the central frequencies {μ₀^∥, μ₀^⊥} of the mother wavelets were chosen to be twice the smallest Fourier scale available, giving {μ₀^∥, μ₀^⊥}={0.023 h Mpc⁻¹, 0.028 h Mpc⁻¹}. This choice allows sampling on the largest possible scales, while still enabling the estimation of statistics on non-periodic fields (see Sect. 4.4). We then chose (J_∥, J_⊥) = (6, 6), corresponding to (Q_∥, Q_⊥) = (0.90, 0.93), which is close to the commonly used dyadic scaling, resulting in a total of 36 wavelets. Finally, we chose the values {σ₀^∥, σ₀^⊥}={0.011 h Mpc⁻¹, 0.015 h Mpc⁻¹} to cover the entire k_∥ − k_⊥ Fourier plane as uniformly as possible. Fig. 2 shows the resulting wavelet set in this cylindrical Fourier space.

Fig. 2. Different regions in the k∥ vs k⊥ plane characterised by wavelets in the 3D wavelet set.

This wavelet set applies not only to EoR lightcones but also to all spectroscopic datasets and to any 2+1D field whose third axis represents a different quantity. In this work, we demonstrate its performance on EoR lightcones with limited redshift evolution (Δz ≃ 1); extending the method to fully evolving lightcones with larger Δz is a direction for future work.

3.2. Scattering transform statistics

Scattering transform (ST) statistics are a family of summary statistics. In this work, we consider a simplified version³ of the scattering covariance statistics introduced in Cheng et al. (2023). Consider a 3D lightcone I. Its scattering covariance statistics are defined as

$$\begin{array}{c}\hfill C_{{\lambda }_1,{\lambda }_2}^{p_1,p_2}=\mathrm{Cov}({[I⊛{\mathrm{\Psi }}_{{\lambda }_1}]}^{p_1},{[I⊛{\mathrm{\Psi }}_{{\lambda }_2}]}^{p_2}),\end{array}$$(9)

where ⊛ denotes a convolution and

$$\begin{array}{c}\hfill {[z]}^p=\{\begin{array}{cc}|z|,\hfill & \mathrm{if}p=0.\hfill \\ z,\hfill & \mathrm{if}p=1.\hfill \end{array}\end{array}$$(10)

The covariance between two stochastic processes X and Y is defined as Cov(X, Y) = ⟨XY^*⟩ − ⟨X⟩⟨Y^*⟩, where ⟨⟩ denotes the expected value and ^* denotes the complex conjugate. In this work, the expected values were estimated using a spatial average.

Using this definition, the three subsets of scattering covariance coefficients are:

$$\begin{array}{ccc}\hfill C^{11}& =\mathrm{Cov}(I⊛{\mathrm{\Psi }}_{{\lambda }_1},I⊛{\mathrm{\Psi }}_{{\lambda }_2}),\hfill & \hfill {\lambda }_1={\lambda }_2,\end{array}$$(11)

$$\begin{array}{ccc}\hfill C^{01}& =\mathrm{Cov}(|I⊛{\mathrm{\Psi }}_{{\lambda }_1}|,I⊛{\mathrm{\Psi }}_{{\lambda }_2}),\hfill & \hfill {\lambda }_1\le {\lambda }_2,\end{array}$$(12)

$$\begin{array}{ccc}\hfill C^{00}& =\mathrm{Cov}(|I⊛{\mathrm{\Psi }}_{{\lambda }_1}|,|I⊛{\mathrm{\Psi }}_{{\lambda }_2}|),\hfill & \hfill {\lambda }_1\le {\lambda }_2.\end{array}$$(13)

In these equations, the condition λ₁ ≤ λ₂ denotes μ_j∥, 2^∥ ≤ μ_j∥, 1^∥ and μ_j⊥, 2^⊥ ≤ μ_j⊥, 1^⊥, as summarised in Fig. 3. This λ-scale conditioning is explained as follows. For C¹¹, the covariance is zero if λ₁ ≠ λ₂. Since the spectral supports of different wavelets do not overlap, I ⊛ Ψ_λ1 and I ⊛ Ψ_λ2 also have distinct spectral supports and therefore exhibit vanishing covariance, as they are linearly decorrelated⁴. For C⁰¹, convolving I ⊛ Ψ_λ1 with the first wavelet concentrates the resulting field into a band of scales around λ₁. Applying the modulus operator, |I ⊛ Ψ_λ1|, demodulates the field by re-centering its support towards lower scales (see, e.g. Zhang & Mallat 2019). To obtain a non-vanishing covariance, the secondary wavelet convolution I ⊛ Ψ_λ2, centred at scales λ₂, must now overlap the envelope defined by |I ⊛ Ψ_λ1|, enforcing λ₁ ≤ λ₂. For C⁰⁰, since the covariance is commutative, |I ⊛ Ψ_λ1| and |I ⊛ Ψ_λ2| play the same role, so the condition λ₁ ≤ λ₂ avoids double counting.

Fig. 3. Schematic of the secondary convolution criterion λ1 ≤ λ2, for the 3D wavelets set in the k∥ vs k⊥ plane. The starting wavelet ψλ1 (black) and the second wavelet ψλ2 probe regions such that μj∥, 2∥ ≤ μj∥, 1∥ and μj⊥, 2⊥ ≤ μj⊥, 1⊥. These wavelets, for the example ψλ1, are shown in the highlighted region (red).

Although the scattering covariances C^p₁, p₂ capture the non-Gaussian properties of a given lightcone, they also depend on the amplitude of its power spectrum. To correct for this dependence, we normalised the scattering covariances by the C¹¹ of a reference field. In this work, this reference field is the target lightcone and the normalised scattering covariance statistic is denoted $C_{\mathrm{ref}}^{11}$. The normalised scattering covariance statistics are defined as follows:

$$\begin{array}{cc}& {\overline{C}}^{00}({\lambda }_1,{\lambda }_2)=\frac{C^{00}({\lambda }_1,{\lambda }_2)}{\sqrt{C_{\mathrm{ref}}^{11}({\lambda }_1)C_{\mathrm{ref}}^{11}({\lambda }_2)}},\hfill \end{array}$$(14)

$$\begin{array}{cc}& {\overline{C}}^{01}({\lambda }_1,{\lambda }_2)=\frac{C^{01}({\lambda }_1,{\lambda }_2)}{\sqrt{C_{\mathrm{ref}}^{11}({\lambda }_1)C_{\mathrm{ref}}^{11}({\lambda }_2)}},\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\overline{C}}^{11}(\lambda )=\frac{C^{11}(\lambda )}{C_{\mathrm{ref}}^{11}(\lambda )}·\hfill \end{array}$$(16)

To further constrain the power spectrum, we used an additional set of C¹¹ coefficients, denoted C^11′. These coefficients were computed according to Eq. (11), with a set of wavelets defined as in Sect. 3.1, but with a different number of scales, J′_α > J_α. These C^11′ coefficients, computed from a larger set of wavelets, used in conjunction with C¹¹, provide an overall improved spectral resolution for power spectrum constraints (Mousset et al. 2024), without increasing the computational complexity for higher-order terms. For these coefficients, we used (J′_∥, J′_⊥) = (15, 15), corresponding to (Q_∥, Q_⊥) = (2.53, 2.60).

We concatenated the various ST coefficients to form the set of ST statistics:

$$\begin{array}{c}\hfill \varphi (I)=\{{\overline{C}}^{11},{\overline{C}}^{11\prime },{\overline{C}}^{00},{\overline{C}}^{01}\}.\end{array}$$(17)

The choice of parameters to define the wavelet sets results in 36 coefficients in C¹¹, 225 coefficients in C^11′, and 441 coefficients in C⁰¹ and C⁰⁰, for a total of 1143 ST statistics.

4. Generative models

4.1. Maximum entropy microcanonical models

We constructed a maximum entropy model of an EoR lightcone from the scattering covariance statistics on the target lightcone, d_t. The new approximate realisations of this lightcone were then sampled using gradient descent (Bruna & Mallat 2018; Allys et al. 2020; Régaldo-Saint Blancard et al. 2023). This gradient descent sampling transforms a white Gaussian distribution into the maximum entropy distribution conditioned on the ST statistics of d_t.

The procedure is as follows. The map on which the gradient descent is performed is denoted u. It is initialised with a white noise realisation, u = u₀, with a mean and standard deviation matching those of the target lightcone. The optimisation is performed by constraining the scattering transform statistic, ϕ(u), of u with the following loss:

$$\begin{array}{c}\hfill \mathcal{L}(u)=\Vert \varphi (d_t)-\varphi (u){\Vert }^2,\end{array}$$(18)

where ∥ ⋅ ∥ is the Euclidean norm and d_t is fixed. This gradient descent is performed until the loss converges. The samples, $\stackrel{\sim }{d}$, from the ST generative model correspond to the map $\stackrel{\sim }{d}=u_f$ obtained at the end of the optimisation, whose ST statistics are close to those of d_t. Different $\stackrel{\sim }{d}$ samples are obtained by initialising the optimisation with different white noise realisations.

4.2. Numerical experiment and statistical validation

To validate the quality of the ST generative models, we compared the statistical properties of the generated lightcones, after performing the inverse quantile transform, with those of the target lightcone, d_t. We generated 30 independent synthesised lightcones, denoted ${\left\{{\stackrel{\sim }{d}}_i\right\}}_{i=1\cdots 30}$ after the inverse quantile transform, each initialised from a different white Gaussian noise realisation. To perform the optimisation, we computed the loss function gradient using automatic differentiation with the L-BFGS-B algorithm (Byrd et al. 1995) (as implemented in PyTorch). We performed the optimisation using 2000 gradient descent iterations, at which point the loss had converged to a plateau after being reduced by approximately ten orders of magnitude. Each optimisation for a single lightcone required two hours on a Nvidia Tesla A100 GPU with 80 GB of memory.

Beyond an initial visual comparison, the statistics considered were the histogram, power spectrum, and Minkowski functionals. In each case, we compared the mean statistics of d_t (shown in blue in comparison figures) and of the 30 ${\stackrel{\sim }{d}}_i$ (shown in green). To assess potential model bias, we compared the difference between these mean statistics and their sample standard deviation, σ_t^oct, between octants of the lightcone⁵, estimated on the target lightcone, d_t. The standard deviation of the St statistics estimator used to construct the generative model is σ_t^oct, since only a volume corresponding to a single octant of the lightcone was used, due to the issue of non-periodic boundary conditions, as explained above. This sampling variance can then simply explain modelling biases of this order.

Additionally, we compared σ_t and σ_s, the sample standard deviation across lightcones, estimated from d_t and the 30 ${\stackrel{\sim }{d}}_i$, respectively. This allowed us to quantify whether the ST generative model reproduces both the mean and sample variance of the statistics. While σ_s can be estimated directly from the 30 ${\stackrel{\sim }{d}}_i$, σ_t has to be estimated from d_t only. To do this, we approximated the octants of d_t as independent, resulting in ${\sigma }_t\simeq {\sigma }_t^{\mathrm{oct}}/\sqrt{8}$. While this is a rough approximation, which will a priori tend to underestimate σ_t, since the octant are not independent, it enables a first comparison between σ_t and σ_s.

To provide additional comparisons, we also considered the distributions of statistics obtained from 30 Gaussian realisations with the same power spectrum as d_t (shown in orange), with their associated lightcone-to-lightcone sample standard deviation, σ_g, and with the distributions obtained from the complete set of lightcones in the Loreli II database (shown in purple).

A list of the different σ values defined in this section is given in Table 1.

4.3. Data conditioning

Prior to the construction of the generative model, the data was conditioned to mitigate several practical issues that could adversely affect the resulting synthesis. This data conditioning ensured that the generative model operated on a well-suited representation of the EoR lightcone, thus improving the fidelity and stability of the subsequent synthesis.

4.3.1. Quantile function

The EoR lightcone is a complex non-Gaussian process, with a highly skewed histogram that is difficult to recover using an ST maximum entropy model. To improve the modelling, we first applied an invertible non-linear transform to the EoR lightcone to obtain a histogram closer to a Gaussian distribution. We then sampled the ST maximum entropy model in this transformed data space, before returning to the original EoR lightcone data space using the inverse non-linear transform.

Such ST modelling based on an invertible non-linear transform has been used in previous studies, for example using a logarithmic pointwise transform before constructing ST generative models (Allys et al. 2020; Régaldo-Saint Blancard et al. 2023; Mousset et al. 2024). This approach is particularly useful for modelling processes with complex histograms, such as those with only non-negative values, including the matter density of large-scale structures in the Universe. The general idea of transforming a process into another one that is better approximated is also used when constructing lognormal models for processes whose logarithm is expected to be much better described by a Gaussian model than the process itself. Lognormal models are used, for example, to describe the thermal emission of Galactic molecular clouds (Miville-Deschênes et al. 2007; Levrier et al. 2018).

In this work, we used a quantile transformation to map the histogram of the target lightcone to a normal distribution (Mark Beasley & Erickson 2009; McCullagh & Tresoldi 2020). To do this, we first estimated the global empirical cumulative distribution function (CDF) $\widehat{F}$ of the field for the entire (δν ∼ 14 MHz) lightcone. We then computed the corresponding percentile p for each data point d(r) in the lightcone,

$$\begin{array}{c}\hfill p(\mathit{r})=\widehat{F}(d(\mathit{r})).\end{array}$$(19)

Next, this percentile p was passed into the inverse of the standard normal CDF (also known as the quantile function) to obtain

$$\begin{array}{c}\hfill \overline{d}(\mathit{r})={\mathrm{\Phi }}^{-1}(p(\mathit{r})).\end{array}$$(20)

This procedure produces a set of $\overline{d}$ values that follow a normal distribution, replacing the original skewed distribution while preserving the ordering of the original data points.

In this work, the quantiles for the transform were estimated from all data points within the target lightcone. The choice of the number of quantiles corresponds to a standard bias-variance trade-off: too few quantiles result in a poor representation of the data distribution, while too many can introduce overfitting and unnecessary computational overhead. We empirically find that 10 000 quantiles provide a good balance. Given that the target lightcone contains approximately 10⁶ pixels, using 10 000 quantiles corresponds to approximately 300 pixels per quantile. This quantile-to-pixel ratio is sufficient to accurately represent the low density tail of the distribution without incurring needless computational cost.

It should be noted that the quantile transform is not a pointwise operator but is instead applied by analysing the entire lightcone simultaneously. Although this transform efficiently models complex histograms, any modelling error in the transformed data space can cause significant instabilities in the inverse transform. In particular, the significant effect of large scales on the histogram, which are difficult to model statistically from a small amount of data, can make it difficult to use such a transform.

4.3.2. Treatment of unconstrained scales

The generative model in this work does not constrain two limiting ranges of scales. The first unconstrained range, denoted κ_low, corresponds to the largest scales, which have Fourier mode wavenumbers smaller than the mother wavelet’s wavenumber, μ₀. The second unconstrained range, denoted κ_high, corresponds to the smallest scales, which have Fourier modes beyond the Nyquist frequency.

Since we first constructed the ST generative model on the quantile-transformed data before applying the inverse quantile transform, these unconstrained scales can significantly affect the final synthesised lightcone because the transform is the non-linear and non-local. To mitigate their impact, we handled them explicitly, as described below:

• The κ_low scales were replaced with the corresponding modes from the target lightcone in the transformed data space. This substitution occurred prior to the gradient descent optimisation, ensuring that these large-scale modes were included in the initial conditions. This was necessary because, if left unconstrained, these scales would influence the global inverse quantile transform. Since these scales are not characterised by the ST statistics, they do not directly influence the generation of the remaining scales during the optimisation process.

• The κ_high scales in the target lightcone were filtered prior to the quantile transform, both when constructing the generative model and when comparing with the synthesised maps. Although the non-linearity of the quantile transform introduces some beyond-Nyquist scales after the histogram transform, their effect on the inverse histogram transform after modelling is negligible.

In principle, our framework can be run without re-inserting the modes with k < k_low from the target lightcone. In this case, these modes remain unconstrained by the optimisation process. The generative model still visually reproduces the target lightcone and accurately matches the scattering statistics on intermediate and small scales. However, the largest scales show biased power, with deviations at k ≲ κ_low comparable to, although generally within, the sample variance of the target lightcone. We also experimented with replacing the modes with Gaussian realisations matched to the target power spectrum. However, because the quantile transform is a global, non-linear mapping, the large-scale modes strongly influence the histogram. Consequently, even modest higher-order mismatches propagate into systematic biases at low k after the inverse transform. A more physical alternative would be to source these large-scale modes from a dedicated low-resolution EoR simulation. However, consistently coupling them to the high-resolution, quantile-transformed small scales is non-trivial, and we leave this as an interesting direction for future work.

4.4. Non-periodic boundary conditions

Simulations of the EoR lightcone have non-periodic boundaries, in particular because of the redshift evolution along the line of sight. However, since we compute convolutions by products in Fourier space, these are inherently periodic. To avoid taking boundary discontinuities into account when estimating the ST statistics, the empirical covariances were computed only over a central smaller cubic volume of the data⁶. This is achieved by cropping the edges of the wavelet convolved fields, so that all pixels over which the average is computed are sufficiently separated from the boundaries. To verify this, the cropped band sizes are matched to the maximum sizes characterised by the wavelets, which are determined by the mother wavelet scale $\frac{1}{{\mu }_0}$. Specifically, we removed a width at each end of each axis equal to half the mother wavelet scale on that axis, yielding an average over a volume of (61, 61, 87) pixels.

Note that while the statistics of the target lightcones, estimated on d_t and used to constrain the generative model, were computed taking into account the non-periodic boundary conditions as described above, the samples generated by this generative model are full lightcones with periodic boundary conditions. To achieve this, their initial conditions were white noise defined over the entire (128, 128, 170) pixels data cube, and their statistics were computed using period convolutions and averages over the entire cube. We performed these syntheses by comparing the two types of ST statistics estimators, with and without periodic boundary conditions, in Eq. (18). As a side effect, the statistics used to construct the generative models were estimated from only about one eighth of the target lightcone volume, leading to a higher sampling variance.

5. Validation of the generative models

5.1. Visual validation

A first test of synthesis quality is to visually compare three realisations, ${\stackrel{\sim }{d}}_i$, of a synthesised lightcone with the target lightcone, d_t. The xy-slice comparison (Fig. 4) shows that these three realisations are visually similar to d_t. The target lightcone contains several distinctive features, such as strongly localised absorbing (negative) regions and ionised (near-zero) regions, which serve as touchstones that the synthesised lightcones successfully capture. Fig. 5 shows that the synthesised lightcones again reproduce the different phases of the field (emission, absorption, and ionised) and the elongated structures characteristic of the redshift axis. However, the generated samples appear to exhibit a greater diversity in the distribution of ionised regions, although a quantitative assessment is difficult from a single d_t lightcone.

Fig. 4. Comparison of a given xy-slice (single frequency channel in the lightcone), for the target EoR lightcone and three synthesised lightcones generated using the 3D scattering transforms with the inverse quantile transform applied. The synthesised lightcones visually reproduce the target lightcone well.

Fig. 5. Same as Fig. 4, but for a given xz-slice (with z representing the redshift domain) for the target EoR lightcone and three synthesised lightcones. As in Fig. 4, the target lightcone is visually well reproduced by the synthesised lightcones.

5.2. Histogram

Fig. 6 compares normalised histograms (i.e. the estimated probability density functions) of the target (blue) and synthesised (green) lightcones. The target lightcone exhibits a non-symmetric, left-skewed tail distribution, a clear, highly non-Gaussian signature, which the synthesised lightcones accurately reproduce over approximately four orders of magnitude, down to δT_b = −7 mK. The bottom plot of Fig. 6 shows that the difference between these histograms, normalised by σ_t^oct, mostly lies between −1 and 1, consistent with the estimation variance of the target ST statistics. The top plot of Fig. 6 shows that σ_s, the lightcone-to-lightcone standard deviation for the synthesised maps, is of the same order as the bin-to-bin variability of the histogram estimated on d_t, indicating that this variance is also well reproduced. The standard deviation, σ_t, is not shown in the top plot, as the sample variability is already apparent due to the fine binning.

Fig. 6. Linear- and log-scale histograms of the target EoR lightcone and a synthesised lightcone after the inverse quantile transform. The histogram of the synthesised lightcone, showing the mean of 30 realisations and their standard deviation (shaded region), recovers the histogram of the target lightcone well, including its skewness.

We used negative brightness temperature values (δT_b) as a proxy for ionised pixels to infer the neutral fraction of the lightcone. The neutral fraction of the target lightcone is 0.617. The mean neutral fraction of the 30 generated lightcones is 0.624 ± 0.001, indicating that the generative model can recover the target lightcone’s neutral fraction to 1% accuracy.

5.3. Power spectrum

Four power spectra were used to compare the target and synthesised lightcones: the 3D isotropic power spectrum; the 2D isotropic power spectrum on the xy or xz planes; and the cylindrically averaged power spectrum. Although the cylindrically averaged power spectrum is sometimes referred to as the 2D power spectrum, here it specifically refers to the power spectrum estimated on actual 2D fields: each 2D xy or xz slice, averaged over the z or y axis, respectively. We refer the reader to Appendix A for more details. Comparisons with Gaussian realisations are not included in this section, as these have the same power spectrum as the target lightcone.

Figure 7a compares the 3D isotropic power spectrum of the target lightcone and that of the synthesised lightcones. The power spectrum of the synthesised lightcones closely follows that of the target lightcone, remaining within ±σ_s on most scales. Although estimating σ_t^oct for the power spectrum is challenging, due to the differing numbers of modes per bin for different data sizes, the averaged 3D power spectrum of the synthesised lightcones reproduces that of d_t within 25% across the range 0.1 and 1 h Mpc⁻¹. Figures 7b and 7c show similar results for the 2D power spectra in the xy− and xz− planes. In these cases, a visible bias appears at the smallest scales, with syntheses reproducing the d_t power spectra to within 15% on average across 0.2−1 h Mpc⁻¹. The shape of the power spectra in Figure 7c results from discretisation effects introduced by the binning scheme in the cylindrical Fourier plane, which imposes a structure on the spectra.

Fig. 7. Comparison between the target and synthesised lightcones after the inverse quantile transform, showing power spectra across different dimensions. (a) Spherically-averaged 3D power spectrum. (b) 2D power spectrum estimated on xy-planes. (c) 2D power spectrum estimated on xz-planes. In all cases the power spectra of the synthesised lightcones reproduce the largest scales well, with increasing deviations at smaller scales.

Fig. 8 compares of the cylindrically averaged power spectra, where the synthesised lightcones are represented by the mean spectrum of the 30 realisations. The right panel shows the ratio of the target lightcone’s cylindrically averaged power spectrum to that of the mean of the synthesised lightcones. Overall, the cylindrical power spectra are well reproduced, but Fig. 8c highlights large ratio values at high k_∥ and k_⊥. These correspond to the bias we see at high k in the power spectrum comparison in Fig. 7, close to Nyquist frequencies. Apart from the bias observed at high k modes, the cylindrical power spectra of the synthesised lightcones remain within $\pm 1.5\stackrel{\sim }{\sigma }$ of the variance across the 30 realisations at all other scales.

Fig. 8. Comparisons of the cylindrically averaged power spectrum of the target lightcone (left) and synthesised lightcones (middle) after the inverse quantile transform. The ratio between the two cylindrically averaged power spectra is also shown (right), with the power spectra of the synthesised lightcones reproducing that of the target lightcone within 25%, on most scales below k∥ = 0.3 h Mpc−1.

The high-k bias in the power spectra of the synthesised lightcones, particularly along the z domain, is not a limitation of the C¹¹ and C^11′ power spectrum constraints. Prior to the application of the inverse quantile transform, power spectra recover well at these scales, with only small bias. However, this bias increases significantly when the transform is subsequently applied. This likely arises because non-Gaussian structures, i.e. the statistical dependence between the different scales, are not properly reproduced at small scales. In this scenario, the highly non-linear inverse histogram transform could introduce such a bias. To address this issue, the C¹¹ and C^11′ power spectrum constraints could be added after this transform, at the expense of additional computational time.

5.4. Minkowski functionals

Finally, we compared the Minkowski functionals (MFs) of the target and the synthesised lightcones, which were computed using the Minkowski functionals function of the QuantImPy Python package (Boelens & Tchelepi 2021). As with the 2D isotropic power spectrum, we computed the 2D MFs on the xy and xz planes, which were then averaged over the z or y axes, respectively. The three MFs are M₀, M₁, and M₂ (Mecke et al. 1994; Schmalzing & Buchert 1997). These 2D MFs were computed by progressively increasing a threshold value and evaluating the properties of the regions above this threshold: total area (M₀), perimeter (M₁), and genus (M₂), which provides insight into the connectivity and complexity of ionised structures.

To address the non-periodic boundary conditions of d_t, the MFs were computed on the xy or xz planes after applying a complete symmetric padding to restore periodicity without significantly altering the underlying data (this results in some sharp edges at the initial boundary position). We also applied this symmetric padding when estimating the MFs of the syntheses, even though they have periodic boundary conditions, to ensure the best possible comparison.

Figures 9 and 10, show the MFs in the xy and xz planes, respectively. Since the lightcone-to-lightcone standard deviations (σ_t, σ_s, and σ_g) are very small, they are multiplied by a factor of 10 in the plots. These plots show that all three MFs of the synthesised lightcones closely resemble those of d_t in both xy and xz planes. This indicates that the ST generative models effectively reproduce the morphological features of the data. However, while the differences between these mean statistics mostly remain within ±σ_t^oct for the xz-MFs, this is not entirely the case for the xy-MFs, for which clear deviations are particularly visible for M₂. We interpret this as the first limitation of the ST generative model, arising from its difficulty in reproducing the sharp, high-contrast, and highly textured circular features present in the target lightcone (see Fig. 4). Overall, for both xy and xz cases, the M₀, M₁, and M₂ statistics are reproduced within 1%, 10%, and 15%, respectively.

Fig. 9. Comparison of the Minkowski functionals applied to each frequency channel of the lightcone between the synthesised lightcones (after the inverse quantile transform) and the target lightcone. The mean Minkowski functionals across the lightcone are shown, with shaded regions indicating their standard deviation. From the bottom panel, M0 and M1 of the synthesised lightcones recover the target lightcone’s Minkowski functionals well, lying within the octant variance, σtoct, and reproducing its morphology. The M2 of the synthesised lightcone follows that of the target lightcone and lies within ±2.5σtoct.

Fig. 10. Same as Fig. 9, but for the xz-plane. The synthesised lightcones reproduce the functionals of the xz-plane, recovering the morphology of the evolving lightcone. These functionals lie within the octant variance, σtoct, of the target lightcone.

Notably, the synthesised data have a lightcone-to-lightcone standard deviation, σ_s, which is a few times smaller than the standard deviation σ_t estimated from d_t. In contrast, Gaussian models exhibit considerably larger variance. This could indicate a second limitation of ST generative models. However, care should be taken when interpreting this result, as estimating σ_t from σ_t^oct could lead to a substantial underestimation of this variance, which could be larger between lightcones than within them.

6. Conclusions

The 21 cm signal from the Epoch of Reionisation (EoR) is observed by radio interferometers as a lightcone that includes a 2D sky plane and a 1D redshift domain that captures the signal’s evolution. These observations are contaminated by foregrounds that are several orders of magnitude stronger than the EoR signal, which the instruments are designed to detect. Recently developed component separation methods rely on maximum entropy generative models based on scattering transforms. In this work, we explored the feasibility of generative models of EoR lightcones built from 3D scattering transforms.

We constructed the 3D scattering transforms from a wavelet set obtained by taking the tensor product between a 1D wavelet set defined in k_z and a 2D wavelet set defined in the k_x − k_y plane. We demonstrated the performance of these generative models by applying them to a single simulated EoR lightcone from the LoReLi II database. We validated the generative model by comparing several realisations of the synthesised lightcone with the original target lightcone, using independent statistical estimators, including histograms, 3D and 2D power spectra, cylindrically averaged power spectra, and Minkowski functionals.

The results demonstrate that the synthesised lightcones accurately reproduce the key non-Gaussian properties and morphology of the target lightcone, including its characteristic skewed histogram, redshift-dependent structures, and morphological complexity, as probed by the Minkowski functionals. These statistics were well reproduced within the expected sample variance estimated from the data.

Beyond the overall validation, several limitations were identified:

• The synthesised lightcones reproduce the target power spectra over a wide range of scales, but clear biases appear at high-k (small scales), particularly along the redshift axis. This indicates limitations of the generative model when combined with the inverse quantile transformation. However, this could be mitigated by additional post-inverse corrections or by refining the treatment of scales not constrained by the scattering transforms.

• The synthesised lightcone Minkowski functionals do not fully reproduce the M₂ of the target lightcone. This suggests that the generative model struggles to produce sharp, high-contrast circular features in the target lightcone, such as ionised regions. More sophisticated scattering transform representation, such as complete scattering spectra statistics (Cheng et al. 2023), could address this, but at an additional computational cost.

• The synthesised lightcones display diverse realisations, but the variance between them may be slightly overestimated compared to the variance estimated over the octants of the target lightcone. More simulated data would be required to adequately estimate this lightcone-to-lightcone variability.

Overall, this study provides the first proof-of-concept for applying 3D scattering transform-based generative models to EoR lightcones. These results pave the way for the development of a statistical component separation method using 3D scattering transform statistics to separate the EoR signal from a realistic forward-modelled lightcone with foregrounds and noise.

Although our tests focus on EoR lightcones with a modest redshift extent, the approach is generally applicable to other spectroscopic datasets and 3D fields in which the third axis represents a different quantity. A natural next step is to apply the method to full lightcones with stronger redshift evolution across the volume.

Acknowledgments

The post-doctoral contract of Ian Hothi was funded by Sorbonne Université in the framework of the Initiative Physique des Infinis (IDEX SUPER) and through the overheads of Edith Falgarone’s MIST ERC. We would like to thank Florent Mertens and Adélie Gorce for their useful discussions and input. The Loreli project was provided with computer and storage resources by GENCI at TGCC thanks to the grants 2022-A0130413759 and 2023-A0150413759 on the supercomputer Joliot Curie’s ROME partition. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. This research made use of astropy, a community-developed core Python package for astronomy (Astropy Collaboration 2018); scipy, a Python-based ecosystem of open-source software for mathematics, science, and engineering (Virtanen et al. 2020) – including numpy (Walt et al. 2011); matplotlib, a Python library for publication quality graphics (Hunter 2007).

References

Appendix A: Power spectrum binning

The dimensionless power spectrum of the 21cm brightness temperature I is defined as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_i^2(k)=\frac{A_k}{2{\pi }_i^2}{\displaystyle \int |\stackrel{\sim }{I}(\mathit{k})·W_i}(\mathit{k}){|}^{2}d\mathit{k},\end{array}$$(A.1)

where k (=k_x, k_y, k_z) is the wavenumber and W_i(k) is the spectral window function that defines the ith bin. We choose the spectral window function to be a top-hat function defined as,

$$\begin{array}{c}\hfill W_i(\mathit{k})=\{\begin{array}{cc}1,\hfill & \mathrm{if}{\mathit{k}}_i\le \mathit{k}<{\mathit{k}}_{i+1}\hfill \\ [1mm]0,\hfill & \mathrm{otherwise}\hfill \end{array}\end{array}$$(A.2)

where k_i defines the lower wavenumber limit of the i^th bin. These window functions are assumed to be isotropic, so W(k) = W(k) with k = |k|. For the 3D power spectrum A_k = k³ and for the 2D power spectrum A_k = k².

Before estimating the power spectrum, it is essential to account for the non-periodic boundary conditions of the target lightcone. We address this by applying a Blackman-Harris window function to the data prior to the Fourier transform. Specifically, for the 3D power spectrum, we apply the 1D Blackman-Harris window function separately along each axis of the 3D lightcone. Similarly, for the 2D power spectrum, this window is applied independently along each axis of the respective 2D field. This windowing procedure is applied consistently to both the target d_t and synthesised ${\stackrel{\sim }{d}}_i$ lightcones; although the ${\stackrel{\sim }{d}}_i$ have periodic boundary conditions, this application ensures consistent treatment of boundary effects.

In addition, we also restrict the power spectra estimation to the range of scales that are well constrained by the generative model, so that the range of k modes that are binned is limited to κ_low ≤ k ≤ κ_high (see Section 4.3.2 for definitions).

The dimensionless cylindrically averaged power spectrum, which bins the k_∥ and k_⊥ domains separately, is defined as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_i^2(k_{\Vert },k_{\perp })=\frac{A_k}{2{\pi }_i^2}{\displaystyle \int |\stackrel{\sim }{I}}({\mathit{k}}_{\Vert },{\mathit{k}}_{\perp })·W_i({\mathit{k}}_{\Vert },{\mathit{k}}_{\perp }){|}^{2}d{\mathit{k}}_{\Vert }d{\mathit{k}}_{\perp },\end{array}$$(A.3)

where A_k = k³, for $k=\sqrt{{{\mathit{k}}_{\Vert }}^2+{{\mathit{k}}_{\perp }}^2}$. We again assume an isotropic top-hat function, W_i(k_∥, k_⊥), for k_∥ = |k_∥| and k_⊥ = |k_⊥|, yielding

$$\begin{array}{c}\hfill W_i(k_{\Vert },k_{\perp })=\{\begin{array}{cc}1,\hfill & \mathrm{if}{k}_{\Vert ,i}\le k_{\Vert }<k_{\Vert ,i+1}\&k_{\perp ,i}\le k_{\perp }<k_{\perp ,i+1}\hfill \\ [1mm]0,\hfill & \mathrm{otherwise}\hfill \end{array}\end{array}$$(A.4)

where k_∥, i defines the lower wavenumber limit of the i^th along the k_∥ domain, and k_⊥, i defines the lower wavenumber limit of the i^th along the k_⊥ domain.

All Tables

All Figures

	Fig. 1. Schematic of the cylindrically-averaged power spectrum of an observed EoR lightcone, with the three distinct regions highlighted.
In the text

	Fig. 2. Different regions in the k∥ vs k⊥ plane characterised by wavelets in the 3D wavelet set.
In the text

	Fig. 3. Schematic of the secondary convolution criterion λ1 ≤ λ2, for the 3D wavelets set in the k∥ vs k⊥ plane. The starting wavelet ψλ1 (black) and the second wavelet ψλ2 probe regions such that μj∥, 2∥ ≤ μj∥, 1∥ and μj⊥, 2⊥ ≤ μj⊥, 1⊥. These wavelets, for the example ψλ1, are shown in the highlighted region (red).
In the text

	Fig. 4. Comparison of a given xy-slice (single frequency channel in the lightcone), for the target EoR lightcone and three synthesised lightcones generated using the 3D scattering transforms with the inverse quantile transform applied. The synthesised lightcones visually reproduce the target lightcone well.
In the text

	Fig. 5. Same as Fig. 4, but for a given xz-slice (with z representing the redshift domain) for the target EoR lightcone and three synthesised lightcones. As in Fig. 4, the target lightcone is visually well reproduced by the synthesised lightcones.
In the text

	Fig. 6. Linear- and log-scale histograms of the target EoR lightcone and a synthesised lightcone after the inverse quantile transform. The histogram of the synthesised lightcone, showing the mean of 30 realisations and their standard deviation (shaded region), recovers the histogram of the target lightcone well, including its skewness.
In the text

Fig. 7. Comparison between the target and synthesised lightcones after the inverse quantile transform, showing power spectra across different dimensions. (a) Spherically-averaged 3D power spectrum. (b) 2D power spectrum estimated on xy-planes. (c) 2D power spectrum estimated on xz-planes. In all cases the power spectra of the synthesised lightcones reproduce the largest scales well, with increasing deviations at smaller scales.

In the text

Fig. 8. Comparisons of the cylindrically averaged power spectrum of the target lightcone (left) and synthesised lightcones (middle) after the inverse quantile transform. The ratio between the two cylindrically averaged power spectra is also shown (right), with the power spectra of the synthesised lightcones reproducing that of the target lightcone within 25%, on most scales below k∥ = 0.3 h Mpc−1.

In the text

Fig. 9. Comparison of the Minkowski functionals applied to each frequency channel of the lightcone between the synthesised lightcones (after the inverse quantile transform) and the target lightcone. The mean Minkowski functionals across the lightcone are shown, with shaded regions indicating their standard deviation. From the bottom panel, M0 and M1 of the synthesised lightcones recover the target lightcone’s Minkowski functionals well, lying within the octant variance, σtoct, and reproducing its morphology. The M2 of the synthesised lightcone follows that of the target lightcone and lies within ±2.5σtoct.

In the text

	Fig. 10. Same as Fig. 9, but for the xz-plane. The synthesised lightcones reproduce the functionals of the xz-plane, recovering the morphology of the evolving lightcone. These functionals lie within the octant variance, σtoct, of the target lightcone.
In the text