The information regarding how the intergalactic medium is reionized by astrophysical sources is contained in the tomographic threedimensional 21 cm images from the epoch of reionization. In Zhao et al. Paper I we demonstrated for the first time that density estimation likelihoodfree inference DELFI can be applied efficiently to perform a Bayesian inference of the reionization parameters from the 21 cm images. Nevertheless the 3D image data needs to be compressed into informative summaries as the input of DELFI by e.g. a trained 3D convolutional neural network CNN as in Paper I DELFI3D CNN. Here in this paper we introduce an alternative data compressor the solid harmonic wavelet scattering transform WST which has a similar yet fixed i.e. no training architecture to CNN but we show that this approach i.e. solid harmonic WST with DELFI outperforms earlier analyses based on 3D 21 cm images using DELFI3D CNN in terms of credible regions of parameters. Realistic effects including thermal noise and residual foreground after removal are also applied to the mock observations from the Square Kilometre Array. We show that under the same inference strategy using DELFI the 21 cm image analysis with solid harmonic WST outperforms the 21 cm power spectrum analysis. This research serves as a proof of concept demonstrating the potential to harness the strengths of WST and simulationbased inference to derive insights from future 21 cm lightcone image data.
['Reionization', 'H I line emission', 'Astrostatistics', 'Bayesian statistics', 'Wavelet analysis', '1383', '690', '1882', '1900', '1918', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
Simulationbased Inference of Reionization Parameters from 3D Tomographic 21 cm Lightcone Images. II. Application of Solid Harmonic Wavelet Scattering Transform
{'Simulation-based Inference of Reionization Parameters from 3D Tomographic 21 cm Light-cone Images - II: Application of Solid Harmonic Wavelet Scattering Transform': "X/i.pc/a.pc/o.pc/s.pc/h.pc/e.pc/n.pc/g.pc Z/h.pc/a.pc/o.pc , 1, 2 Y/i.pc M/a.pc/o.pc , 1 S/h.pc/i.pc/f.pc/a.pc/n.pc Z/u.pc/o.pc , 3, 4 /a.pc/n.pc/d.pc B/e.pc/n.pc/j.pc/a.pc/m.pc/i.pc/n.pc D. W/a.pc/n.pc/d.pc/e.pc/l.pc/t.pc 2, 5, 6  \n1 Department of Astronomy, Tsinghua University, Beĳing 100084, China  \n2 Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris (IAP), 98 bis bd Arago, 75014 Paris, France  \n3 National Astronomical Observatories, Chinese Academy of Sciences, Beĳing 100101, China  \n4 Key Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beĳing, 100101, P. R. China  \n5 Sorbonne Université, Institut Lagrange de Paris (ILP), 98 bis bd Arago, 75014 Paris, France  \n6 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA", 'ABSTRACT': "The information regarding how the intergalactic medium is reionized by astrophysical sources is contained in the tomographic three-dimensional 21 cm images from the epoch of reionization. In Zhao et al. (2022a) ('Paper I'), we demonstrated for the first time that density estimation likelihood-free inference (DELFI) can be applied efficiently to perform a Bayesian inference of the reionization parameters from the 21 cm images. Nevertheless, the 3D image data needs to be compressed into informative summaries as the input of DELFI by, e.g., a trained 3D convolutional neural network (CNN) as in Paper I ( DELFI-3D CNN ). Here in this paper, we introduce an alternative data compressor, the solid harmonic wavelet scattering transform (WST), which has a similar, yet fixed (i.e. no training), architecture to CNN, but we show that this approach (i.e. solid harmonic WST with DELFI) outperforms earlier analyses based on 3D 21 cm images using DELFI-3D CNN in terms of credible regions of parameters. Realistic effects, including thermal noise and residual foreground after removal, are also applied to the mock observations from the Square Kilometre Array (SKA). We show that under the same inference strategy using DELFI, the 21 cm image analysis with solid harmonic WST outperforms the 21 cm power spectrum analysis. This research serves as a proof of concept, demonstrating the potential to harness the strengths of WST and simulation-based inference to derive insights from future 21 cm light-cone image data.  \nKeywords: Reionization (1383), H I line emission (690), Astrostatistics (1882), Bayesian statistics (1900), Wavelet analysis (1918)", '1. INTRODUCTION': "The intensity mapping of the 21 cm line associated with the spin-flip transition of H I atoms is a promising probe of the epoch of reionization (EoR; Furlanetto et al. 2006). It contains information regarding when and how the intergalactic medium (IGM) was heated and reionized by the first luminous objects. Upper limits of the 21 cm power spectrum from the EoR (Paciga et al. 2013; Parsons et al. 2014; Jacobs et al. 2015; Mertens et al. 2020; Trott et al. 2020; Yoshiura et al. 2021; Abdurashidova et al. 2022) have been placed by ongoing interferometric array experiments, including the Precision Array for Probing the Epoch of Reionization (PAPER; Parsons et al. 2010), the Murchison Widefield Array (MWA;  \nCorresponding author: Xiaosheng Zhao, Yi Mao  \nzhaoxs18@mails.tsinghua.edu.cn (XZ), ymao@tsinghua.edu.cn (YM)  \nTingay et al. 2013), the LOw Frequency ARray (LOFAR; van Haarlem et al. 2013), and the Giant Metrewave Radio Telescope (GMRT; Intema et al. 2017). In the foreseeable future, the first measurements of the 21 cm power spectrum from the EoR will be very likely achieved by upcoming array experiments including the Hydrogen Epoch of Reionization Array (HERA; DeBoer et al. 2017) and the Square Kilometre Array (SKA; Mellema et al. 2013). Furthermore, the SKA will have the exciting promise of mapping the three-dimensional (3D) tomographic light-cone images of the 21 cm brightness temperature from the EoR with high sensitivity.  \nThe21cmsignalis non-Gaussian due to reionization patchiness. Therefore, the 3D light-cone images of the 21 cm signal contain more information than the power spectrum statistic, which is unlike the traditionally well-studied case of Gaussian signals in the cosmic microwave background (CMB). As such, it is of key importance for 21 cm observers to understand how to optimally extract the information in the 3D 21 cm images.  \nThe 3D image data needs to be compressed into informative summaries inevitably because it is technically very challenging to process the high-dimensional image data directly. For this purpose, several works have proposed to apply the convolutional neural networks (CNNs) to compress the 2D 21 cm image slices (Gillet et al. 2019) or 3D 21 cm light-cone images (Zhao et al. 2022a, hereafter referred to as 'Paper I'; Prelogović et al. 2022; Neutsch et al. 2022) into data summaries. However, the practical applications of CNN-based methods are generically computationally expensive both in generating a large volume of simulation data that are needed for training the networks and in the process of training and optimization itself. Even with so much 'engineering', the fine-tuned networks may still be sub-optimal (see, e.g. Paper I; Prelogović et al. 2022).  \nTo mitigate these problems, it has been proposed to inject the inductive bias into CNNs with scattering transform (Mallat 2012a; Allys et al. 2019; Cheng et al. 2020; Pedersen et al. 2022) and utilize the scattering transform to construct scattering or wavelet networks (Gauthier et al. 2021; Pedersen et al. 2022). The scattering transform employs the filters that have well-behaved mathematical structures, e.g. the Morlet filters (Mallat 2012a; Trott 2016), and under its unique definition, exploits the modulus nonlinearities and hierarchical structures - similar to the multi-layers in the CNN - which allows it to extract information across multi-scales. Compared with CNNs, the scattering transform has fixed filter parameters, and therefore do not need to be trained, which is a significant advantage against the CNNs. Recently, the scattering transform has been extended to 3D applications. For example, the harmonic-related wavelets are introduced to infer the molecular properties (Eickenberg et al. 2017, 2018), cosmological parameters (Saydjari et al. 2021; Valogiannis & Dvorkin 2022; Chung 2022; Eickenberg et al. 2022), and the CMB B-mode (Jeffrey et al. 2022). Specifically, Eickenberg et al. (2022) employs the first-order wavelet-based features and shows the advantage of harmonic wavelets against the isotropic and oriented ones. In the context of extracting the information from the 3D 21 cm light-cone (i.e. 'lightcuboids'), in this paper, we apply the solid harmonic wavelet scattering transform (WST; Eickenberg et al. 2017, 2018) to compress the 3D image data. The solid harmonic WST injects the inductive bias into 3D CNNs with both 3D solid harmonic wavelets and the scattering transform which outputs multipleorder wavelet-based features.  \nTraditionally, the posterior distributions for the parameters of the reionization model (hereafter referred to as 'reionization parameters') can be inferred from the measurements of statistical observables of the 21 cm signal - i.e. informative summaries from the data scientific point of view - with the Monte Carlo Markov Chain (MCMC) analysis that uses an explicit (e.g. Gaussian) likelihood approximation (see the  \n21CMMC code; Greig & Mesinger 2015, 2017, 2018). However, the inference using the explicit likelihood approximation in 21 cm analysis may find itself either biased if the assumption in likelihood function is not exact or the covariance matrices in the likelihood neglect the off-diagonal elements between different wavenumbers and different redshifts, or computationally too expensive if all elements of covariance matrices are properly accounted for (Zhao et al. 2022b, hereafter referred to as Z22b; Prelogović & Mesinger 2023). To address this issue of intractable likelihood, the so-called 'implicit likelihood inference' (ILI; Alsing et al. 2018, 2019; Papamakarios 2019; Cranmer et al. 2020; Tejero-Cantero et al. 2020), aka 'simulation-based inference' (SBI) or 'likelihoodfree inference' (LFI), has been recently proposed to 'learn' the density of the likelihood or posterior directly from data, using advanced methods in deep learning, e.g. conditional masked autoregressive flows (CMAFs; Papamakarios et al. 2017) which is a variant of normalizing flows (Papamakarios et al. 2021).  \nIn Paper I, we introduced the density estimation likelihoodfree inference (DELFI) to the 21 cm analysis for the first time and performed the posterior inference of reionization parameters from the 3D tomographic 21 cm light-cone images. The 3D CNNs were adopted therein to compress the 3D 21 cm images into informative summaries. In this paper, we improve the data compressor of 3D 21 cm images and replace the 3D CNNs in Paper I with the solid harmonic WST /one.sup , but otherwise still perform the Bayesian inference of reionization parameters in the framework of DELFI. This new approach (i.e. solid harmonic WST with DELFI) is dubbed ' 3D ScatterNet ' herein. We will compare the inference results using 3D ScatterNet and those using DELFI-3D CNN , in order to demonstrate the improvement of the data compressor. In addition, we will compare the inference results from the 21 cm images using 3D ScatterNet and those from the 21 cm power spectrum analysis using 21cmDELFI-PS (Z22b), both in the framework of DELFI, in order to demonstrate whether solid harmonic WST can extract more information from the 21 cm images than the power spectrum statistic. Similar to the analysis in Z22b, realistic effects, including thermal noise and residual foreground after removal with the singular value decomposition (SVD; Stewart 1993; Wolz et al. 2015; Masui et al. 2013), are applied to the mock observations from the", 'SKA.': 'The rest of this paper is organized as follows. The 3D ScatterNet methodology is introduced in Section 2. We describe the simulations and the application of realistic effects in Section 3, present the inference results in Section 4, and  \nmakeconcludingremarksinSection5. Sometechnical details are left to Appendix A (on the effect of angular frequency), Appendix B (on the effect of smoothing scales), Appendix C (on the dependence of light-cone volumes), Appendix D (on the network setting and sample size), and Appendix E (on the performance of DELFI-3DCNN with SKA noise). Some of our results were previously summarized by us in a conference paper (Zhao et al. 2023b). Our code is publicly available on GitHub /github .', '2.1. Solid Harmonic Wavelet Scattering Transform': "We briefly summarize the solid harmonic WST in this subsection, following Eickenberg et al. (2017, 2018).  \nSolid harmonic WST convolves the original fields with a cascade of solid harmonic wavelets, performs non-linear moduli on these convolved fields, and integrates over the coordinate space. The solid harmonic wavelet is defined as the solid harmonic multiplying a Gaussian  \n𝜓 𝑚 ℓ ( u ) = 1 ( √ 2 𝜋 ) 3 𝑒 -| 𝑢 | 2 / 2 | u | ℓ 𝑌 𝑚 ℓ GLYPH<18> u | u | GLYPH<19> , (1)  \nwhere | u | ℓ 𝑌 𝑚 ℓ GLYPH<16> u | u | GLYPH<17> is the solid harmonic, evaluated at the coordinate u . The Gaussian part serves as localizing the wavelet support around zero. Here we omit the further normalization factors for simplicity. In order to capture features of the field at multiple scales, the mother solid harmonic wavelet in Equation (1) is dilated at the scale /two.sup 2 𝑗 , i.e.  \n𝜓 𝑚 𝑗,ℓ ( u ) = 2 -3 𝑗 𝜓 𝑚 ℓ GLYPH<0> 2 -𝑗 u GLYPH<1> . (2)  \nThe Euclidean norm as the modulus operator, also dubbed ' first-order modulus coefficient' herein, is defined by  \n𝑈 [ 𝑗 , ℓ ] d ( u ) = ℓ ∑︁ 𝑚 = -ℓ GLYPH<12> GLYPH<12> GLYPH<12> d ∗ 𝜓 𝑚 𝑗,ℓ ( u ) GLYPH<12> GLYPH<12> GLYPH<12> 2 ! 1 / 2 , (3)  \nwhere the field d is convolved (denoted by ' ∗ ' in Equation 3) by the dilated solid harmonic wavelets at the scale 𝑗 with the angular frequency band ℓ . The additional rotational phase subspace information represented by 𝑚 is aggregated to produce coefficients that are covariant to rotation. The translational covariance is also guaranteed by the convolution operation.  \nThe first-order solid harmonic wavelet scattering coefficient is given by  \n𝑆 1 [ d ; 𝑗 , ℓ, 𝑞 ] = ∫ R 3 | 𝑈 [ 𝑗 , ℓ ] d ( u )| 𝑞 𝑑 3 u , (4)  \nwhere the modulus is raised by the power 𝑞 , which results in the coefficients that are sensitive to the amplitude of the field, in that the small (large) value of 𝑞 gives more weight to the small (large) non-zero values of the integral. Note that the point-wise transformation by the power 𝑞 does not change the covariance property. By integrating these coefficients over the position u , we get the coefficients that are invariant to both translation and rotation.  \nIn order to capture the information across multiple scales, the first-order modulus coefficient 𝑈 [ 𝑗 , ℓ ] d ( u ) is convolved with another wavelet at a different scale 𝑗 ' with 𝑗 ' > 𝑗 but with the same angular frequency band ℓ , i.e. the second-order modulus coefficient is defined as  \n𝑈 [ 𝑗 , 𝑗 ' , ℓ ] d ( u ) = ℓ ∑︁ 𝑚 = -ℓ GLYPH<12> GLYPH<12> GLYPH<12> 𝑈 [ 𝑗 , ℓ ] d ∗ 𝜓 𝑚 𝑗 ' ,ℓ ( u ) GLYPH<12> GLYPH<12> GLYPH<12> 2 ! 1 / 2 , 𝑗 < 𝑗 ' . (5)  \nThe covariance property is also maintained.  \nThe second-order solid harmonic wavelet scattering coefficient is defined by integration over the coordinate space, similar to Equation (4),  \n𝑆 2 [ d ; 𝑗 , 𝑗 ' , ℓ, 𝑞 ] = ∫ R 3 | 𝑈 [ 𝑗 , 𝑗 ' , ℓ ] d ( u )| 𝑞 𝑑 3 u , (6)  \nwhich is also invariant to both translation and rotation.  \nIn principle, successive similar operations can be applied to define the higher-order solid harmonic wavelet scattering coefficients /three.sup . In this paper, we neglect the information encoded in the higher-order scattering coefficients and only consider the zeroth-, first- and second-order scattering coefficients, following Allys et al. (2019); Cheng et al. (2020). In the same spirit as in Equations (4) and (6), the zeroth-order scattering coefficient is defined as the sum of all pixel values raised by the power 𝑞 , i.e.  \n𝑆 0 [ d ; 𝑞 ] = ∫ R 3 d ( u ) 𝑞 𝑑 3 u . (7)  \nBasically, while the zeroth-order scattering coefficients highlight the amplitude of the field d , the first-order scattering coefficients separate the information of the field into different scales by gathering the specific geometric feature in the original field at the scale 𝑗 . Also, the second-order scattering coefficients represent the information of nonlinear mode mixing between different scales 𝑗 and 𝑗 ' . These nonlinear operations obviously extract a richer set of information than the power spectrum analysis that only characterizes the features at separate Fourier modes. In addition, not only does the WST capture the local spatial information with the  \nFigure 1. Visualization of solid harmonic WST. Shown is one slice of the original field of simulated cosmological 21 cm signals along the angular direction (top row), the first-order modulus coefficients 𝑈 [ 𝑗 , ℓ ] d ( u ) (middle row) and the second-order modulus coefficients 𝑈 [ 𝑗 , 𝑗 ' , ℓ ] d ( u ) (bottom row), at some representative wavelet scales 𝑗 and/or 𝑗 ' and at various angular frequencies ℓ . The colormap of all fields is the same as shown for the original field.  \n<!-- image -->  \nFigure 2. Thefirst-order (' S 1 ') and second-order (' S 2 ') scattering coefficients of the simulated cosmological 21 cm signals and their dependence on the reionization parameters. The scattering coefficients are averaged over 0 ≤ ℓ ≤ ℓ max and evaluated at 𝑞 = 1. We vary two reionization parameters separately, namely, log 10 ( 𝑇 vir ) (top) and log 10 𝜁 (bottom) - see their definitions in Section 3.1 - with the fiducial model in green lines.  \n<!-- image -->  \nlocalized wavelets, but also it saves some large-scale information because of the long tail of the wavelets, similar to the Morlet transform (Trott 2016). Furthermore, the scattering coefficients are naturally invariant to translation and rotation. Also, they are Lipschitz continuous to deformations (Mallat 2012b), meaning that they are approximately proportional to small deformations of the original field.  \nSolid harmonic WST is very analogous to 3D CNN in three aspects as follows. Firstly, the modulus operators work similarly to the nonlinear functions in the latter. Secondly, the integration of the modulus coefficients over the coordinate space is essentially the pooling operation. Thirdly, the orders in scattering coefficients are similar to the layers in the CNN. However, solid harmonic WST has fixed kernels in wavelets, so unlike the CNN, it does not need to be trained in order to output the scattering coefficients. Also, large kernels are usually avoided in the CNN because otherwise, this would require a huge set of training parameters, but large kernels can be applied to solid harmonic WST efficiently.  \nFigure 1 presents the visualization of selected modulus coefficients. Each coefficient reflects some features of the original field. For example, the first-order modulus coefficients with ( 𝑗 = 0, ℓ = 1) highlight the boundaries of H II bubbles, while the second-order modulus coefficients with ( 𝑗 = 0, 𝑗 ' = 1, ℓ = 1) show more complex structures. Also, as ℓ increases (from the left to right in Figure 1), the modulus coefficients encode the information of small structures of the original field. Note that wavelets with ℓ = 0 are Gaussian, while wavelets with nonzero ℓ can encode rich structures such as filaments.  \nTo compute the scattering coefficients, we choose the maximum scale 𝑗 max = 5 and maximum angular frequency ℓ max = 6, and select the value of the modulus power 𝑞 = 0 . 5, 1, or 2 for the first and second-order scattering coefficients. The half-width parameter (i.e. standard deviation in the Gaussian factor) of the mother solid harmonic wavelets is set to be unity, in which case the maximum scale 𝑗 max = 5, which satisfies the criterion 2 × ( half width parameter ) × 2 𝑗 max ≤ the number of simulation cells on each side (66 in this work, as shown in Section 3.1). For a given ℓ (with 0 ≤ ℓ ≤ ℓ max) and 𝑞 , the number of the first order coefficients is 𝑗 max + 1 = 6 since 𝑗 takes the value of 0 ≤ 𝑗 ≤ 𝑗 max and the number of the second order coefficients is 𝑗 max ( 𝑗 max + 1 )/ 2 = 15 since 𝑗 and 𝑗 ' take the value of 𝑗 = 0, 1 ≤ 𝑗 ' ≤ 𝑗 max; 𝑗 = 1, 2 ≤ 𝑗 ' ≤ 𝑗 max; etc. To keep the dimension of coefficients reasonably low, in this paper, we simply average the information over different ℓ for a given 𝑞 . (We will discuss the effect of angular frequencies ℓ on the parameter inference in Appendix A.) Therefore, the first and second-order coefficients have a total of 63 components for three values of power 𝑞 = 0 . 5, 1, and 2.  \nThezeroth-order coefficient is ill-defined if 𝑞 = 0 . 5because the field value at a single pixel can be negative. Also, the zeroth-order coefficient is zero if 𝑞 = 1 because it is simply the global sum of the field in this case which is zero for the 21 cm signal by definition (see Section 3.1). Therefore, instead of 𝑞 = 0 . 5 or 1, we choose (arbitrarily) three higher modulus power 𝑞 = 2, 3, and 4 in this paper for the zerothorder coefficients.  \nAltogether, the final concatenated coefficient vector has a dimension of 66, which can be adjusted if necessary. Furthermore, following Allys et al. (2019), we take the logarithms of those coefficients with base 2 as our new coefficients. (If a component is negative, we take the logarithm of the absolute value but keep the sign of the component.) The advantage of adopting these new coefficients is that they behave linearly with 𝑗 and 𝑗 ' . Hereafter the phrase 'scattering coefficients' refers to their logarithms.  \nWeshowthe representative values of scattering coefficients in Figure 2. We vary the reionization parameters (see their definitions in Section 3.1), and find that the coefficients are generically suppressed (enhanced) as log 10 𝑇 vir (log10 𝜁 ) increases. This implies that the scattering coefficients are sensitive to the reionization parameters.", '2.2. Simulation-based inference with CMAFs': 'Figure 3 shows the workflow of 3D ScatterNet . The scattering coefficients extracted by the solid harmonic WST serve as the input summaries of the CMAFs /four.sup , which is a variant of the neural density estimators (NDEs) that perform the implicit likelihood inference. We refer interested readers to Paper I for an in-depth description of the DELFI and the CMAFs. The detailed settings of NDEs are given in Appendix D. In this paper, we set two neural layers of a single transform and 50 neurons per layer for CMAFs. We also use the ensembles of CMAFs to improve the performance. The number of transforms and the configuration of ensembles are chosen based on the performance of posterior validation.', '2.3. Validation of the posterior': "Asthe final step of inference, we validate both marginalized and joint posteriors. In the posterior validation, hypothesis tests are made to check if the posteriors from CMAFs are self-consistent statistically. Note that the posterior validation is sometimes referred to as 'posterior calibration'. However, this step is not to calibrate the trained networks, but simply an indication of whether the network complexity and training data are adequate to learn the conditional density accurately in a statistical manner. We follow Appendix A of Paper I and recap the validation statistics in this subsection.  \nFigure 3. The workflow of 3D ScatterNet . (Left, fixed part) the data compression with solid harmonic WST - the 3D light-cone images d which are simulated with the parameters 𝜽 are compressed by a cascade of scattering transforms (each containing the convolution with solid harmonic wavelets, harmonic modulus, and integration operation) to form the summaries t which include the zeroth-, first-, and second-order scattering coefficients { S 0 , S 1 , S 2 } , respectively. (Right, trainable part) The CMAFs are used to learn the summary density conditional on the parameters (i.e. likelihood). With the likelihood, the posterior can be inferred using Bayes' Theorem at the data summary t 0 .  \n<!-- image -->  \nFor the marginalized posterior, the probability integral transform (PIT; Gneiting et al. 2007; Mucesh et al. 2021) is defined as  \nPIT ( ˜ 𝜃 ) ≡ ∫ ˜ 𝜃 -∞ 𝑓 ( 𝜃 ) d 𝜃 , (8)  \ni.e. the cumulative distribution function (CDF) of the inferred marginal distribution 𝑓 ( 𝜃 ) at the true value ˜ 𝜃 . If the inferred posteriors are accurate, then the PIT should be uniformly distributed.  \nFor the joint posteriors, there are two statistics. The first one is the copula probability integral transformation (copPIT; Ziegel & Gneiting 2014; Mucesh et al. 2021) which is defined as  \ncopPIT ≡ Pr { 𝐻 ( 𝜽 ) ≤ 𝐻 ( ˜ 𝜽 )} . (9)  \nHere 'Pr' represents the probability and 𝐻 ( 𝜽 ) is the CDF of the inferred joint distribution. The copPIT is the multivariate extension of PIT.  \nThesecondstatistic is the highest probability density (HPD; Harrison et al. 2015) which is defined as  \nHPD ( ˜ 𝜽 ) ≡ ∫ 𝑓 ( 𝜽 ) ≥ 𝑓 ( ˜ 𝜽 ) 𝑓 ( 𝜽 ) d 𝑛 𝜽 . (10)  \nThe HPD describes the plausibility of ˜ 𝜽 under the distribution 𝑓 ( 𝜽 ) . A small value of HPD indicates high plausibility. Similar to the PIT, the copPIT and HPD should be also uniformly distributed if the posteriors are accurate.  \nIn order to check the uniformity, we adopt two metrics, the Kolmogorov-Smirnov (KS; Kolmogorov 1992) test and Cramér-von Mises (CvM; Anderson 1962) test, which focus on different aspects of distribution. While the KS test is  \nsensitive to the median, the CvM test captures the tails of a distribution. If the 𝑝 -value from a test is larger than the preset value of significance level, typically 0.01 or 0.05 (Željko Ivezić et al. 2014), then the null hypothesis that these statistics follow a uniform distribution is accepted. Throughout this paper, we adopt the significance level of 0.01 and report the results only if the 𝑝 -value is larger than 0.01 unless stated otherwise.", '3.1. Cosmic 21 cm Signal': "The 21 cm brightness temperature relative to the CMB temperature at position x can be written (Furlanetto et al. 2006) as  \n𝑇 21 ( x , 𝑧 ) = ˜ 𝑇 21 ( 𝑧 ) 𝑥 HI ( x ) [ 1 + 𝛿 ( x )] GLYPH<18> 1 -𝑇 CMB 𝑇 𝑆 GLYPH<19> , (11)  \nwhere ˜ 𝑇 21 ( 𝑧 ) = 27 √︁ [( 1 + 𝑧 )/ 10 ] ( 0 . 15 / Ω m ℎ 2 )( Ω b ℎ 2 / 0 . 023 ) in units of mK, 𝑥 HI ( x ) is the neutral fraction, and 𝛿 ( x ) is the matter overdensity, at position x . We assume the cold dark matter can be traced by the baryon perturbation on large scales, so 𝛿 𝜌 H = 𝛿 . In this paper, we focus on the limit where spin temperature 𝑇 𝑆 ≫ 𝑇 CMB, likely valid soon after reionization begins, though this assumption is strongly model-dependent. As such, we can neglect the dependence on spin temperature. Also, as a demonstration of concept, we ignore the effect of peculiar velocity; such an effect can be readily incorporated in forward simulations by the algorithm introduced by Mao et al. (2012).  \nFrequency[MHz]  \nFigure 4. Visualization of the mock 21 cm observation in a slice along the LOS. From top to bottom: (1) the simulated cosmic 21 cm images without realistic effect applied except that the k ⊥ = 0 mode is removed from each 2D slice perpendicular to the LOS (i.e. 'pure signal'); (2) same as the pure signal but with total noise from SKA observation (i.e. 'SKA noise'); (3) same as the pure signal but with total noise from SKA observation and the foreground contamination; (4) same as the pure signal but with total noise from SKA observation and residual foreground after the foreground is removed with the SVD technique (i.e. 'SKA noise + residual foreground'). The map of pure signal retains the original grid size in simulations, but the rest of the three maps are made by convolving the original maps with a Gaussian filter in the angular direction and a top-hat filter in the LOS direction, with the widths of the filters corresponding to the size of the 1-km baseline.  \n<!-- image -->  \nIn this paper, we use the publicly available code 21cmFAST /five.sup (Mesinger & Furlanetto 2007; Mesinger et al. 2011), which can be used to perform semi-numerical simulations of reionization, as the simulator to generate the datasets. Our simulations were performed on a cubic box of 100 comoving Mpc on each side, with 66 3 grid cells. Following the interpolation approach in Paper I, the snapshots at nine different redshifts of the same simulation box (i.e., with the same initial condition) are interpolated to construct a light-cone 21 cm data cube within the comoving distance of a simulation box along the line of sight (LOS). We concatenate 10 such light-cone boxes, each simulated with different initial conditions in density fields but with the same reionization parameters, together to form a full light-cone datacube /six.sup of the size 100 × 100 × 1000 comoving Mpc 3 (or 66 × 66 × 660 grid cells) in the redshift range 7 . 51 ≤ 𝑧 ≤ 11 . 67. To mimic the observations from radio interferometers, we subtract from the light-cone field the mean of the 2D slice for each 2D slice perpendicular to the LOS, because radio interferometers cannot measure the mode with k ⊥ = 0.  \nWeparametrize our reionization model as follows, and refer interested readers to Paper I for a detailed explanation of their physical meanings.  \n- (1) 𝜁 , the ionizing efficiency, which is a combination of several parameters related to ionizing photons. In our paper, we vary 𝜁 as 10 ≤ 𝜁 ≤ 250.  \n(2) 𝑇 vir, the minimum virial temperature of halos that host ionizing sources. In our paper, we vary this parameter as 4 ≤ log 10 ( 𝑇 vir / K ) ≤ 6.  \nCosmological parameters are fixed in this paper as ( Ω Λ , Ω m, Ω b, ns, 𝜎 8, ℎ ) = (0.692, 0.308, 0.0484, 0.968, 0.815, 0.678) (Planck Collaboration et al. 2016).", '3.2. Thermal noise and residual foreground': 'To apply the realistic effects to the 21 cm signal, we first compute the 𝑢𝑣 coverage map at each of a total of 660 frequency channels. The 𝑢𝑣 coverage map can be used to suppress the thermal noise and calculate the telescope response to the cosmic 21 cm signal and foreground. We then generate the thermal noise with a 1000-hour observation. For the thermal noise estimation in this paper, we employ the Tools21cm /seven.sup (Giri et al. 2020) code to simulate the expected thermal noise in the 3D 21 cm light-cone images as observed with the SKA1-Low. We list our basic assumptions of the  \nTable 1. Specifications for the SKA  \nSKA configuration in Table 1. When synthesizing the total signal, we applied the same telescope response in the calculation of thermal noise to both the cosmic 21 cm signal component and the radio foreground component. The signal of a pixel in the 𝑢𝑣 space is kept unchanged for non-zero 𝑢𝑣 coverage but otherwise set to zero if the pixel is outside the 𝑢𝑣 coverage. In order to suppress the noise, we smooth the map with a relatively small baseline of 1 km /eight.sup , which roughly corresponds to the size of the core area of SKA1Low. Specifically, we convolve the map with a Gaussian filter in the angular direction and a top-hat filter in the LOS direction. The FWHM of the Gaussian filter corresponds to the size of the 1-km baseline and the width of the top-hat filter also matches that FWHM.  \nTo mock the foreground contamination, we employ the pygsm /nine.sup package that is based on the GSM-building model (Zheng et al. 2017), which interpolates the sky maps with 29 sky map observations using an improved principal component analysis (PCA) method. We then use these foreground maps to interpolate our images on the grid for each frequency slice. In order to prevent overfitting, we assign a random patch of sky (except for the North and South poles) to model each foreground image. Finally, we employ the SVD technique for foreground removal. Specifically, we remove the largest six singular value modes.  \nFor visualization purposes, we show the mock 21 cm observation with realistic effects in Figure 4. When the cosmic 21 cm signal is applied with the thermal noise of SKA and/or the residual foreground, the features at small scales are washed out but those at the very large scales are still retained.', '4. RESULTS': "For each approach of data compression, we first perform the inference for two representative mock observations - the 'Faint Galaxies Model' and the 'Bright Galaxies Model', defined in Tables 2 and 3 (see the 'True value' therein), respectively, following the convention of Greig & Mesinger  \n(2017) (see, also, Paper I and Z22b). They are selected as two illustrative examples with extreme parameter values which are nevertheless fine-tuned for these models to have similar global reionization histories. We further perform the posterior validation of the trained NDEs using another set of 300 independent samples that are randomly drawn from the allowed region in the parameter space in which the mean neutral fraction satisfies 0 . 08 ≤ ¯ 𝑥 HI ≤ 0 . 81 at 𝑧 = 7 . 1, corresponding to the 2 𝜎 credible region as constrained by the IGM damping wing of ULASJ1120+0641 (Greig et al. 2017).", '4.1. 3D ScatterNet vs 3D CNN': "We compare the results using 3D ScatterNet and using DELFI-3D CNN in this subsection. For simplicity, here we only consider the case of cosmological 3D 21 cm images without realistic effect applied except that the k ⊥ = 0 mode is removed from each 2D slice perpendicular to the LOS (i.e. 'pure signal'). Figure 5 shows the credible regions (1 𝜎 and 2 𝜎 , or 0.68 and 0.95 levels) for the 'Faint Galaxies Model' and the 'Bright Galaxies Model', with quantitative comparisons in terms of medians and 1 𝜎 (16th and 84th percentile) errors presented in Tables 2 and 3.  \nThe results from DELFI-3DCNN were taken from Paper I, in which the ILI was performed with data compression made by a trained 3D CNN. Paper I used a set of 9,000 samples for training and validation of the 3D CNN and another set of 9,000 samples for training and validation of the density estimators. Nevertheless, if the training sample size was doubled, the results for the 'Faint Galaxies Model' and the 'Bright Galaxies Model' were not improved, which implies that the inference accuracy might be limited by some intrinsic properties in experimental choices (including the network architecture, characteristics in training datasets, and other hyper-parameter choices) for training the 3D CNN. However, when we replace the data compression from 3D CNN to solid harmonic WST, we find that the inference results are significantly improved in terms of the location and size of credible regions in the posterior distributions for the reionization parameters, as shown in Figure 5 and Tables 2 and 3. Also, the degeneracy of these two parameters is clearly revealed in the credible region for 3D ScatterNet , which indicates that this degeneracy is intrinsic in the theoretical modeling.  \nNow that the solid harmonic WST represents a better approach to data compression than 3D CNN, we will focus on the 3D ScatterNet in the remainder of this paper.", '4.2. 3D ScatterNet vs 21cmDELFI-PS': 'Z22b shows that for the power spectrum analysis, the DELFI approach (i.e. 21cmDELFI-PS ) outperforms the standard MCMC analysis. In this subsection, therefore, we compare the results using 3D ScatterNet and 21cmDELFI-PS , both under the same inference strategy using DELFI. On this  \nTable 2. Bayesian Inference with 21cmDELFI-PS , DELFI-3D CNN and 3D ScatterNet for the \'Faint Galaxies Model\'  \nNote. Here, \'pure signal\' refers to the mock observations of cosmological 3D 21 cm images (i.e. without thermal noise or foreground contamination, but the k ⊥ = 0 mode is removed from each 2D slice perpendicular to the LOS); \'SKA noise\' refers to the mock SKA observations of the 3D 21 cm images with total noise (with the contributions from thermal noise and sample variance errors) yet without foreground contamination; \'SKA noise + residual foreground\') refers to the mock SKA observations of the 3D 21 cm images with total noise and residual foreground after the foreground is removed with the SVD technique.  \nTable 3. Same as Table 2 but for the \'Bright Galaxies Model\'Table 4. Recovery Performance by the 21cmDELFI-PS and 3D ScatterNet  \n𝑎 The coefficient of determination 𝑅 2 is computed for the medians (in the logarithmic scale) of the inferred posteriors from 300 testing samples.  \n𝑏 The fractional error 𝜖 refers to the relative error of the deduced parameters in the linear scale. Here we present the 68% credible interval of the probability density distribution of 𝜖 .  \nTable 5. The 𝑝 -values for the Null Hypotheses That These Statistics Are of a Uniform Distribution with 300 Testing Samples  \nsame footing, we can investigate whether data compression using solid harmonic WST extracts more information from the 3D 21 cm light-cone images than the power spectrum statistic.  \n4.2.1. Pure Signal  \nWefirst consider the \'pure signal\' case for the comparison, with the credible regions for the \'Faint Galaxies Model\' and the \'Bright Galaxies Model\' shown in Figure 5 and quantitative results listed in Tables 2 and 3. For the \'Faint Galaxies Model\', the systematic shift (i.e. relative errors of the predicted medians with respect to the true values) and the 1 𝜎 statistical errors are 0 . 04% ± 0 . 13% (0% ± 0 . 23%) for  \n<!-- image -->  \nFigure 5. The posteriors estimated from 3D cosmic 21 cm light-cone images (i.e., without thermal noise or residual foreground) by three different approaches of data compression - power spectrum using 21cmDELFI-PS (green), 3D CNN using DELFI-3D CNN (blue) and solid harmonic WST using 3D ScatterNet (red) - for two mock observations: the \'Faint Galaxies Model\' (left) and the \'Bright Galaxies Model\' (right). We show the 1 𝜎 (dark) and 2 𝜎 (light) credible regions. The dashed lines indicate the true parameter values.  \n<!-- image -->  \nlog 10 ( 𝑇 vir ) with 3D ScatterNet ( 21cmDELFI-PS ), respectively /one.sup/zero.sup , and 0 . 14% + 0 . 47% -0 . 41% (0 . 27% ± 0 . 74%) for log10 𝜁 with 3D ScatterNet ( 21cmDELFI-PS ), respectively. The estimated statistical errors using 3D ScatterNet are about 1 . 7 times smaller than using 21cmDELFI-PS . The similar results hold generically for the \'Bright Galaxies Model\', too.  \nNext, we test the trained NDEs on 300 samples. Table 4 shows the coefficient of determination, 𝑅 2 = 1 -"GLYPH<0> 𝑦 pred -𝑦 true GLYPH<1> 2 / "GLYPH<0> 𝑦 true -𝑦 true GLYPH<1> 2 , where 𝑦 true and 𝑦 pred are the true value and the predicted median in a sample for the parameter 𝑦 (e.g. log 10 ( 𝑇 vir ) and log10 𝜁 in this paper), respectively, and the summation is over all testing samples. 𝑦 true is the average of the true value over all testing samples. A score of 𝑅 2 close to unity indicates an overall good inference performance of this parameter. For the case of \'pure signal\', both 3D ScatterNet and 21cmDELFI-PS give very high 𝑅 2 score, but the 3D ScatterNet slightly outperforms the 21cmDELFI-PS .  \nTable 4 also presents the 68% credible interval of the probability density distribution of 𝜖 , the fractional errors of the  \ndeduced parameter 𝑥 in the linear scale (e.g. 𝑇 vir and 𝜁 in this paper). The fractional error refers to that of the predicted median with respect to the true value, i.e. 𝜖 = ( 𝑥 pred -𝑥 true )/ 𝑥 true. For the case of \'pure signal\', the typical fractional error (represented by the length of the interval) of 3D ScatterNet is about 1.50 times smaller than that of 21cmDELFI-PS for 𝑇 vir, and 1.75 times smaller for 𝜁 .  \nLastly, we perform the posterior validation. Table 5 shows the 𝑝 -value for some hypothesis tests using 300 testing samples. Note that all 𝑝 -values are larger than 0.01 and most are larger than 0.5, which implies that our results are at least reliable with a significance of 0.01. Also, the 𝑝 -values of 3D ScatterNet are generically larger than those of 21cmDELFIPS .', '4.2.2. SKA Noise and Residual Foreground': "We now consider the case of 'SKA noise' which refers to the mock SKA observations of the 3D 21 cm images with total noise (with the contributions from thermal noise and sample variance errors) yet without foreground contamination. Since the images after smoothing with 1-km baseline lose small-scale information, we discard the components of scattering coefficients with 𝑗 = 0 and hence the final scattering coefficients have the dimension of 48. For the power spectrum, we also discard the largest 𝑘 -modes in a single box and generate the final power spectrum vector with a dimension of 120.  \n<!-- image -->  \nFigure 6. Same as Figure 5, but the estimations are made from mock observations of the 3D 21 cm light-cone images from SKA, which includes the total noise (with the contributions from thermal noise and sample variance errors) yet without foreground contamination. Here we only compare two different approaches of data compression - power spectrum using 21cmDELFI-PS (green) and solid harmonic WST using 3D ScatterNet (red).  \n<!-- image -->  \nFor the case of 'SKA noise', Figure 6 shows the credible regions for the 'Faint Galaxies Model' and the 'Bright Galaxies Model', with quantitative comparisons in terms of medians and 1 𝜎 errors presented in Tables 2 and 3. The inference results of the 'SKA noise' case have larger uncertainties than those of 'pure signal', which is reasonable, given the noise. The statistical errors of 3D ScatterNet are smaller than those of 21cmDELFI-PS . This improvement is also verified by Table 4: the 𝑅 2 values of 3D ScatterNet are larger than those of 21cmDELFI-PS , and the typical fractional error of the former is about 1 . 49 (1 . 50) times smaller for 𝑇 vir ( 𝜁 ) than that of the latter.  \nLastly, we consider the case of 'SKA noise + residual foreground' which refers to the mock SKA observations of the 3D 21 cm images with total noise and residual foreground after the foreground is removed with the SVD technique. For the same reason of smoothing as in the 'SKA noise' case, we reduce the dimension of the scattering coefficients to be 48 by discarding the components with 𝑗 = 0. For the power spectrum, in order to train reliable NDEs, we have to further discard more large𝑘 modes with the upper limit 𝑘 max = 0 . 45 Mpc -1 and the final vector of power spectrum has the dimension of 70. In this case, the information for parameter inference is only from large-scale modes because the images are smoothed with a rather coarse resolution.  \nFor the case of 'SKA noise + residual foreground', Figure 7 shows the inference results for two mock observations, with the medians and 1 𝜎 errors listed in Tables 2 and 3. As expected, given more uncertainties due to residual foreground, the posteriors have larger errors in this case than those of 'SKA noise'. The statistical errors of 3D ScatterNet are on average smaller than those of 21cmDELFI-PS . This improvement is also verified by Table 4: the 𝑅 2 values of 3D ScatterNet are higher than those of 21cmDELFI-PS , and the typical fractional error of the former is about 1 . 09 (1 . 15) times smaller for 𝑇 vir ( 𝜁 ) than that of the latter.  \nRegarding the posterior validation, all 𝑝 -values are larger than 0.01 both for the 'SKA noise' case and the case of 'SKA noise + residual foreground', which implies that our results are at least reliable with a significance of 0.01. Also, the 𝑝 -values of 3D ScatterNet are generically larger than those of 21cmDELFI-PS .  \nIn sum, for all cases of different assumptions in noise and residual foreground, 3D ScatterNet outperforms the 21cmDELFI-PS . Our results demonstrate that the solid harmonic WST can extract the information from the 3D 21 cm light-cone images in a more effective way than the power spectrum statistic itself. These results are consistent with the previous findings that the 𝑛 th-order scattering transform captures the information up to the 2 𝑛 -point function (Mallat 2012a; Cheng & Ménard 2021). In addition, Sui et al.  \n(2023) introduced the concept of mutual information to evaluate the comparisons between different data summaries by giving them scores and found that the solid harmonic WST not only outperforms the power spectrum but also surpasses the combined information including both power spectrum and bispectrum, in the context of extracting the information from the 21 cm datacube.", '4.3. Discussions': "In this subsection, we attempt to provide some insights regarding why the 3D CNN (trained with the particular set of experimental choices in Paper I) performs less effectively than the solid harmonic WST (and power spectrum) in compressing the 3D 21 cm light-cone images, and the limitations of the 3D CNN configurations in Paper I.  \nThe solid harmonic WST is analogous to a 3D CNN, but unlike the latter, it essentially gives a fixed kernel without the training process. Furthermore, the resulting scattering coefficients are naturally invariant to translations and rotations and are Lipschitz continuous to deformations (Mallat 2012b). These properties are particularly useful for 3D 21 cm light-cone images because the data summaries have the translational invariance in the angular direction and the rotational invariance along the LOS /one.sup/one.sup , and can be stable under slight deformations due to variations in the initial density fields. These intrinsic invariance properties can enhance conditional density learning of the solid harmonic WST.  \nIn principle, thanks to the similarity between solid harmonic WSTand3DCNN,adeepand/orwide3DCNNshould perform at least as well as the WST in extracting information, as long as three conditions are all met - training data is adequate, network is sufficiently expressive, and training process is successful. Our results imply that the 3D CNN configurations in Paper I might have limitations in these three aspects.  \nThe first limitation is insufficient training data. A dataset that samples the complete parameter-data distribution is required to train a robust neural network. For training a 2D CNN (Gillet et al. 2019), a dataset comprising 9,000 training samples might suffice, but the same number of training samples is likely not adequate for 3D CNN /one.sup/two.sup . While we find that the results for the 'Faint Galaxies Model' and the 'Bright Galaxies Model' were not improved if the training sample size was doubled, a thorough scaling test on the number of training samples is necessary to make a robust diagnosis.  \nThe second limitation is the insufficient complexity of neural network architectures. While the solid harmonic WST  \nexploits its intrinsic invariance properties, the standard CNNs cannot because they are not invariant to rotations due to inherent architecture. Thus the architecture in CNNs needs to be fine-tuned during the training. Alternatively, a likely more effective approach is to integrate additional inductive biases into the architecture (Kauderer-Abrams 2017; Semih Kayhan &van Gemert 2020; Weiler et al. 2018a,b). In these variants of CNNs, the invariance properties can be enhanced, which might optimize the network architecture. Also, in Neutsch et al. (2022), an elongated kernel is designed to detect the patterns evolving in the redshift direction, which reduces the number of trainable parameters /one.sup/three.sup .  \nThe third limitation is the sub-optimal network training process. Network optimization involves refined practices for neural network initialization and learning rate determination (DeZoort & Hanin 2023). Some generic optimization strategies (Neutsch et al. 2022; Akiba et al. 2019) may be useful as well.  \nNote that the 3D CNN trained in Paper I performs even less effectively than the power spectrum analysis using 21cmDELFI-PS . This might be explained by the limitations of 3D CNN discussed above. While the 3D CNN configurations in Paper I might be improved in these three aspects, we leave a thorough exploration along these lines to future work.  \nThe comparison involving 3D CNNs in this work has been limited to the case of pure signal. This is because the comparison when instrumental effects are applied would require additional fine-tuning efforts for the 3D CNN to fully assess its potential, which is not the focus of this work. Appendix E provides some insight into such a comparison by including the SKA noise in the application of DELFI-3DCNN , nevertheless without additional fine-tuning, i.e. using the same hyperparameters for training the 3D CNN as detailed in Paper I.", '5. SUMMARY': "In this paper, we introduce the solid harmonic WST for compressing the 3D image data and generating meaningful low-dimensional summaries. We apply this technique to the data compression of 3D tomographic 21 cm light-cone images and use the resulting scattering coefficients as the input summaries of the DELFI. With DELFI, we perform the Bayesian inference of the reionization parameters where the likelihood is implicitly defined by the forward simulations. This new technique, dubbed 3D ScatterNet (i.e. solid harmonic WST with DELFI), recovers accurate posterior distributions for the reionization parameters.  \nWe compare the inference results with two different approaches of data compression, the solid harmonic WST and  \n<!-- image -->  \nFigure 7. Same as Figure 6, but with mock observations of the 3D 21 cm light-cone images from SKA, which includes the total noise and residual foreground after the foreground is removed using SVD.  \n<!-- image -->  \nthe 3D CNN, for two representative mock observations in the 'pure signal' case, and demonstrate that the 3D ScatterNet outperforms the DELFI-3D CNN significantly. Our results imply that the solid harmonic WST extracts the information from the 3D 21 cm light-cone images in a more informative manner than a 3D CNN given reasonable finetuning. We also offer insights into enhancing the 3D CNN regarding the datasets, network architecture, and training process. Moreover, the solid harmonic WST has fixed kernels in wavelets, which means that it does not need to be trained in order to output the data summaries. This highlights its robustness and efficiency, which is another advantage of the solid harmonic WST over the 3D CNN.  \nWe then make another comparison of the inference results between data compression using the solid harmonic WST ( 3D ScatterNet ) and using the power spectrum statistic ( 21cmDELFI-PS ), both under the same inference strategy using DELFI. The comparisons were made for three cases with different assumptions on the noise and residual foreground 'pure signal', 'SKA noise' and 'SKA noise + residual foreground'. We find that the 3D ScatterNet outperforms the 21cmDELFI-PS in all cases. This implies that the summaries compressed with the solid harmonic WST contain more (i.e. non-Gaussian) information from the 3D 21 cm light-cone images than the power spectrum analysis.  \nOur results demonstrate that combining a WST (3D solid harmonic WST in this paper) with the simulation-based inference will be a promising tool for the scientific interpretation  \nof future 21 cm light-cone image observation data. Based on the findings of this paper, there is room for improvement with regard to the design of summary statistics from solid harmonic WST or new WST. For example, instead of treating the light-cone as a whole to build a statistic, a variant approach might be to treat discrete boxes separately and concatenate the scattering coefficients of each box together as the data summaries. In this way, the correlations between different stages of reionization can be exploited by the CMAFs. Nevertheless, this will increase the dimension of data summaries and therefore incur higher computational costs, so the trade-off between accuracy and efficiency in the solid harmonic WST is yet to be further explored. Also, further cross-correlation between scales resulting from re-scaling the phase information is likely to be exploited in a 3D extension of the wavelet phase harmonics (WPH; Mallat et al. 2019; Allys et al. 2020) for the 2D case. On the other hand, when scaling the similar analysis to a high-dimensional parameter and/or data space, the diffusion model (Luo 2022; Legin et al. 2023; Zhao et al. 2023a) is a promising alternative to the CMAFs, or other normalizing flows, to perform the simulation-based inference because of its simplified training objective. We leave the exploration of these directions to future works.  \nNote. -In the final stage of preparing this manuscript, Greig et al. (2022); Eickenberg et al. (2022) were posted on arXiv. These papers applied the scattering transform (or wavelet moments, varied band-limited first-order scattering coefficients defined in the Fourier space, in Eickenberg et al.  \n2022) for parameter inference, but both used the Fisher matrix formalism. In comparison, we apply the scattering transform for parameter inference using simulation-based inference in this paper. In addition, Greig et al. (2022) focuses on the reionization parameter estimation using the 2D 21 cm images; Eickenberg et al. (2022) focuses on the cosmological parameter estimation using the 3D cosmological density fields. In comparison, our paper has a distinct focus, i.e. on the reionization parameter estimation using the 3D 21 cm light-cone images.", 'ACKNOWLEDGEMENTS': 'This work is supported by the National SKA Program of China (grant No. 2020SKA0110401), NSFC (grant No. 11821303), and the National Key R&D Program of China (grant No. 2018YFA0404502). BDW acknowledges support from the Simons Foundation. We thank Sihao Cheng, Paulo Montero-Camacho, and Bohua Li for their useful discussions  \nand help. We acknowledge the Tsinghua Astrophysics HighPerformance Computing platform at Tsinghua University for providing computational and data storage resources that have contributed to the research results reported within this paper.  \nSoftware: 21CMMC (Greig & Mesinger 2015, 2017, 2018), 21cmFAST (Mesinger & Furlanetto 2007; Mesinger et al. 2011), pydelfi (Alsing et al. 2019), TensorFlow (Abadi et al. 2016), GetDist (Lewis 2019), NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), scikit-learn (Pedregosa et al. 2011), Python2 (Van Rossum & Drake Jr 1995), Python3 (Van Rossum & Drake 2009), Kymatio (Andreux et al. 2018), Tools21cm (Giri et al. 2020), pygsm (de Oliveira-Costa et al. 2008; Zheng et al. 2017), Mayavi (Ramachandran & Varoquaux 2011), galpro (Mucesh et al. 2021), seaborn (Waskom 2021), Astropy (Astropy Collaboration et al. 2013, 2018), healpy (Zonca et al. 2019), HEALPix (Górski et al. 2005).', 'A. THE EFFECT OF ANGULAR FREQUENCY ℓ': "In the main text of this paper, the scattering coefficients are averaged over 0 ≤ ℓ ≤ ℓ max with ℓ max = 6. In this section, we explore the effect of the angular frequency ℓ on the parameter inference.  \nFigure 8 shows the scattering coefficients evaluated at given values of ℓ . While the first-order scattering coefficients show tilting when varying ℓ , the second-order scattering coefficients decrease in overall amplitude as ℓ increases. Also, for ℓ = 0 (i.e. Gaussian wavelets), the scattering coefficients are flat at the scale of 𝑗 ' .  \nWeconstruct two sets of experiments of posterior inference over 300 testing samples and show the coefficient of determination 𝑅 2 of the predicted medians in Figure 9. In the first set of experiments (' ℓ single'), the scattering coefficients are evaluated at each single value of ℓ . Wefindthatusingthe ℓ = 0 information has the smallest values of 𝑅 2 . When using the information at larger ℓ , the 𝑅 2 value increases and then slightly decreases. This implies that the solid harmonic wavelet has better performance than the Gaussian wavelet (i.e. the case with ℓ = 0). In the second set of experiments (' ℓ average'), the scattering coefficients are averaged over 0 ≤ ℓ ' ≤ ℓ for a given value of ℓ . Comparing these two scenarios, we find that the combined information (' ℓ average') gets a higher value of 𝑅 2 generically than the single information (' ℓ single'). When increasing the value of ℓ , in the case of ' ℓ average', the 𝑅 2 value increases and then slightly decreases, with the peak at ℓ = 6. For this reason, in the main text of this paper, we choose the ' ℓ average' scenario with ℓ max = 6.  \nTable 6. Recovery Performance by the 3D ScatterNet for Smoothing with the Baseline of 1-km and 2-km Respectively", 'B. THE EFFECT OF SMOOTHING SCALES': 'In the main text of this paper, the images are smoothed with the scale corresponding to the size of the 1 km baseline. In this section, we test the effect of smoothing scales on parameter inference. We compute the values of 𝑅 2 and typical fractional error 𝜖 for the predicted medians of 300 testing samples. Table 6 shows that smoothing with the 1-km baseline has a larger 𝑅 2 value and smaller fractional error than smoothing with the 2-km baseline. This is because the smoothed images in the former have higher S/N although in coarser resolution. This result is the reason why we chose the size of the 1-km baseline for smoothing in this paper.', 'C. THE DEPENDENCE OF LIGHT-CONE VOLUMES': "In the main text of this paper, the full light-cone images, which concatenate ten light-cone boxes, are exploited for parameter inference. In this section, we investigate whether the information from one of these boxes dominates over the  \nFigure 8. The first-order (' S 1 ') and second-order (' S 2 ') scattering coefficients of the simulated cosmological 21 cm signals, evaluated at given values of ℓ and at 𝑞 = 1.  \n<!-- image -->  \nTable 7. Recovery Performance for Information in Discrete BoxesNote. The parameter inference is made for the case of 'pure signal'. Here, 'First', 'Third' and 'Fifth' refer to the scenario where only the information within the first, third, and fifth box (counting from the low-redshift, corresponding to the central redshift of 7 . 67, 8 . 36, and 9 . 11) is exploited for parameter inference, respectively. 'Full' refers to the case of exploiting the full light-cone images (i.e. concatenating all ten boxes), which is the case in the main text of this paper. † We do not show these results for the first box because the hypothesis test with the HPD values failed, likely due to the lack of training samples.  \nothers. (If this were the case, then observations could be optimized targeting at the 'sweet-spot' redshift.)  \nFor this purpose, we perform the parameter inference from discrete boxes, both using the power spectrum and using scattering coefficients as summaries. Figure 10 shows the posteriors for two mock observations and Table 7 shows the results of recovery performance for the predicted medians of 300 testing samples. For the 3D ScatterNet , the box at the low-redshift results in better inference performance in terms of larger 𝑅 2 , smaller typical fractional error, and smaller credible region for mock observations than that at the high-redshift. However, no single box appears to dominate the information, because the full light-cone images always yield the optimal result that is significantly better than any of the single boxes. These results hold similarly for the 21cmDELFI-PS . We conclude that there is no optimal redshift. As such, the full light-cone images, which gather the complete information, provide the parameter inference with the best performance.", 'D. THE NDE SETTING AND SAMPLE SIZE': "In this section, we give the details of the networks as follows. We choose the CMAFs as the NDEs throughout of this paper. In all architectures, we set two neural layers of a single transform, represented by the masked autoencoders for density estimation (MADE; Germain et al. 2015), with 50 neurons per layer. We also use the ensembles of NDEs. The final posterior is the stacked one from individual posteriors with weights according to their training errors. We fine-tune the training sample size and the CMAFs architecture for every method and experiment based on the outcomes of posterior validation (calibration). The guiding principle behind this fine-tuning is that data (summaries) with greater uncertainties might require more intricate CMAFs, such as increased transformations within a single CMAF and the increased number of CMAFs within the ensemble. Meanwhile, data (summaries) of higher dimensions might necessitate a more substantial training sample due to the vast feature space.  \nThe procedure for scattering coefficients is as follows.  \nFigure 9. The coefficient of determination 𝑅 2 as a function of the angular frequency ℓ . The phrase ' ℓ single' refers to the scenario where the scattering coefficients are evaluated at a given single value of ℓ , while the phrase ' ℓ average' refers to the scenario where the scattering coefficients are averaged over 0 ≤ ℓ ' ≤ ℓ for a given value of ℓ . The parameter after this phrase (i.e. 𝑇 vir or 𝜁 ) indicates the parameter of the 𝑅 2 calculation.  \n<!-- image -->  \n- (1) The scattering coefficients for the 'pure signal' case with k ⊥ = 0 mode removed (see Section 4.1). We use 18,000 light-cone images and an ensemble of 8 NDEs ( 5 , 6 , 7 , 8 ) ∗ 2, which means that we have two NDE blocks, each having the NDE with the number of transformations 5, 6, 7, and 8, respectively. Hereafter we use the same terminology. Notice that the final performance can be further enhanced with double the size of the light-cone images. To fairly compare with 21cmDELFI-PS , we use 18,000 images for reports in the paper.\n- (2) The coefficients from pure discrete boxes (see Appendix C). We use 27,000 samples. The ensembles for the first, third, and fifth boxes: ( 5 , 6 , 7 , 8 ) ∗ 3, ( 5 , 6 , 7 , 8 ) ∗ 3, ( 5 , 6 , 7 , 8 ) ∗ 2, respectively.\n- (3) The coefficients with the information of averaged ℓ from pure boxes (see Appendix A). We use 36,000 samples. The ensembles from ℓ = 1 to ℓ = 8: ( 5 ) ∗ 4, ( 6 , 7 , 8 , 9 ) ∗ 3, ( 6 , 7 , 8 , 9 ) ∗ 2, ( 5 , 6 , 7 , 8 ) ∗ 2, ( 6 , 7 , 8 , 9 ) ∗ 3, ( 5 , 6 , 7 , 8 ) ∗ 2, ( 5 ) ∗ 4, ( 5 , 6 , 7 , 8 ) ∗ 2, respectively.\n- (4) The coefficients with the information of single ℓ from pure boxes (see Appendix A). We use 36,000 samples.  \nThe ensembles from ℓ = 0 to ℓ = 8: ( 5 ) ∗ 4, ( 5 ) ∗ 4, ( 5 , 6 , 7 , 8 ) ∗ 2, ( 5 ) ∗ 4, ( 5 , 6 , 7 , 8 ) ∗ 2, ( 5 , 6 , 7 , 8 ) ∗ 2, ( 5 , 6 , 7 , 8 )∗ 2, ( 5 , 6 , 7 , 8 )∗ 2, ( 5 , 6 , 7 , 8 )∗ 2, respectively.  \n- (5) The coefficients from the signal with SKA noise. We use 27,000 samples (see Section 4.2.2). The ensembles for the signal with 1-km smoothing: ( 5 , 6 , 7 , 8 ) ∗ 2, and 2-km smoothing: ( 5 , 6 , 7 , 8 ) ∗ 3, respectively.\n- (6) The coefficients from the signal with SKA noise and residual foreground (see Section 4.2.2). We use 36,000 samples. The ensembles for the signal with 1-km smoothing: ( 20 ) ∗ 2, and 2-km smoothing: ( 30 ) ∗ 2, respectively.  \nThe procedure for the power spectrum is as follows.  \n- (1) The power spectrum from pure signal with k ⊥ = 0 mode removed (see Section 4.1). We use 18,000 samples and the ensembles ( 5 ) ∗ 4.\n- (2) The power spectrum from pure discrete boxes (see Appendix C). For the first box, we use 32,596 samples and the ensembles ( 6 , 7 , 8 , 9 ) ∗ 3; for the third box, we use 34,044 samples and the ensembles ( 5 , 6 , 7 , 8 ) ∗ 3; and for the fifth box, we use 27,000 samples and the ensembles ( 5 , 6 , 7 , 8 ) ∗ 2.\n- (3) The power spectrum from the signal with SKA noise (see Section 4.2.2). We use 27,000 samples and the ensembles ( 5 ) ∗ 4.\n- (4) The power spectrum from the signal with SKA noise and residual foreground (see Section 4.2.2). We use 36,000 samples and the ensembles ( 20 ) ∗ 3.", 'E. PERFORMANCE OF DELFI-3DCNN WITH SKA NOISE': 'In this section, we extend the comparison of 3D CNN using DELFI-3D CNN , power spectrum using 21cmDELFI-PS and solid harmonic WST using 3D ScatterNet to the mock observations of the 3D 21 cm light-cone images from SKA, which includes the total noise yet without foreground contamination. Ideally, this would require additional fine-tuning efforts for the 3D CNN to fully assess its potential. Here we present the results of DELFI-3D CNN without additional fine-tuning, i.e. using the same hyper-parameters for training the 3D CNN and NDEs as detailed in Paper I. Figure 11 shows that in this realistic scenario, a trained 3D CNN still performs sub-optimal compared to 21cmDELFI-PS and 3D ScatterNet .', 'REFERENCES': "Abadi, M., Barham, P., Chen, J., et al. 2016, in Proceedings of the 12th USENIX Conference on Operating Systems Design and Implementation, OSDI'16 (Berkeley, CA: USENIX  \nAbdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 925, 221, doi: 10.3847/1538-4357/ac1c78  \nFigure 10. The posteriors estimated from 3D cosmic 21 cm light-cone images (i.e. 'pure signal') for the scenario where only the information within the first (red), third (green), fifth (blue) box and the full (black) light-cone images is exploited for parameter inference, respectively. Shown are the results using 21cmDELFI-PS (top) and using 3D ScatterNet (bottom) for the 'Faint Galaxies Model' (left) and the 'Bright Galaxies Model' (right), respectively.  \n<!-- image -->  \n<!-- image -->  \nFigure 11. Same as Figure 6 but including the results from DELFI-3DCNN .  \n<!-- image -->  \nAkiba, T., Sano, S., Yanase, T., Ohta, T., & Koyama, M. 2019, arXiv e-prints, arXiv:1907.10902,  \ndoi: 10.48550/arXiv.1907.10902  \n- Allys, E., Levrier, F., Zhang, S., et al. 2019, A&A, 629, A115, doi: 10.1051/0004-6361/201834975\n- Allys, E., Marchand, T., Cardoso, J. F., et al. 2020, PhRvD, 102, 103506, doi: 10.1103/PhysRevD.102.103506\n- Alsing, J., Charnock, T., Feeney, S., & Wandelt, B. 2019, MNRAS, 488, 4440, doi: 10.1093/mnras/stz1960\n- Alsing, J., Wandelt, B., & Feeney, S. 2018, MNRAS, 477, 2874, doi: 10.1093/mnras/sty819\n- Anderson, T. W. 1962, The Annals of Mathematical Statistics, 33, 1148 , doi: 10.1214/aoms/1177704477\n- Andreux, M., Angles, T., Exarchakis, G., et al. 2018, arXiv e-prints, arXiv:1812.11214. https://arxiv.org/abs/1812.11214 Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,  \nA&A, 558, A33, doi: 10.1051/0004-6361/201322068  \n- Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f\n- Bianco, M., Giri, S. K., Iliev, I. T., & Mellema, G. 2021, Monthly Notices of the Royal Astronomical Society, 505, 3982,  \ndoi: 10.1093/mnras/stab1518  \n- Cheng, S., & Ménard, B. 2021, arXiv e-prints, arXiv:2112.01288. https://arxiv.org/abs/2112.01288\n- Cheng, S., Ting, Y.-S., Ménard, B., & Bruna, J. 2020, MNRAS, 499, 5902, doi: 10.1093/mnras/staa3165  \nChung, D. T. 2022, MNRAS, 517, 1625,  \ndoi: 10.1093/mnras/stac2662  \n- Cranmer, K., Brehmer, J., & Louppe, G. 2020, Proceedings of the National Academy of Sciences, 117, 30055, doi: 10.1073/pnas.1912789117\n- de Oliveira-Costa, A., Tegmark, M., Gaensler, B. M., et al. 2008, MNRAS, 388, 247, doi: 10.1111/j.1365-2966.2008.13376.x\n- DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, Publications of the Astronomical Society of the Pacific, 129, 045001. http://stacks.iop.org/1538-3873/129/i=974/a=045001\n- DeZoort, G., & Hanin, B. 2023, arXiv e-prints, arXiv:2306.11668, doi: 10.48550/arXiv.2306.11668\n- Eickenberg, M., Exarchakis, G., Hirn, M., & Mallat, S. 2017, in Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS'17 (Red Hook, NY, USA: Curran Associates Inc.), 6543-6552\n- Eickenberg, M., Exarchakis, G., Hirn, M., Mallat, S., & Thiry, L. 2018, The Journal of chemical physics, 148, 241732\n- Eickenberg, M., Allys, E., Moradinezhad Dizgah, A., et al. 2022, arXiv e-prints, arXiv:2204.07646,  \ndoi: 10.48550/arXiv.2204.07646  \nFurlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Physics Reports, 433, 181, doi: 10.1016/j.physrep.2006.08.002  \nZhao, X., Zuo, S., & Mao, Y. 2023b, arXiv e-prints, arXiv:2307.09530, doi: 10.48550/arXiv.2307.09530 Zheng, H., Tegmark, M., Dillon, J. S., et al. 2017, MNRAS, 464, 3486, doi: 10.1093/mnras/stw2525 Ziegel, J. F., & Gneiting, T. 2014, Electronic Journal of Statistics, 8, 2619 , doi: 10.1214/14-EJS964  \nZonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298  \nŽeljko Ivezić, Connolly, A. J., VanderPlas, J. T., & Gray, A. 2014, Statistics, Data Mining, and Machine Learning in Astronomy: A Practical Python Guide for the Analysis of Survey Data (Princeton University Press), doi: doi:10.1515/9781400848911"}