Scattering spectra models for physics

Abstract Physicists routinely need probabilistic models for a number of tasks such as parameter inference or the generation of new realizations of a field. Establishing such models for highly non-Gaussian fields is a challenge, especially when the number of samples is limited. In this paper, we introduce scattering spectra models for stationary fields and we show that they provide accurate and robust statistical descriptions of a wide range of fields encountered in physics. These models are based on covariances of scattering coefficients, i.e. wavelet decomposition of a field coupled with a pointwise modulus. After introducing useful dimension reductions taking advantage of the regularity of a field under rotation and scaling, we validate these models on various multiscale physical fields and demonstrate that they reproduce standard statistics, including spatial moments up to fourth order. The scattering spectra provide us with a low-dimensional structured representation that captures key properties encountered in a wide range of physical fields. These generic models can be used for data exploration, classification, parameter inference, symmetry detection, and component separation.


INTRODUCTION
An outstanding problem in statistics is to estimate the probability distribution () of high dimensional data  from few or even one observed sample.In physics, establishing probabilistic models of stochastic fields is also ubiquitous, from the study of condensed matter to the Universe itself.Indeed, even if physical systems can generally be described by a set of differential equations, it is usually not possible to fully characterize their solutions.Complex physical fields, described here as non-Gaussian random processes , may indeed include intermittent phenomena as well as coherent geometric structures such as vortices or filaments.Having realistic probabilistic models of such fields however allows for considerable applications, for instance to accurately characterize and compare non-linear processes, or to separate different sources and solve inverse problems.Unfortunately, no generic probabilistic model is available to describe complex physical fields such as turbulence or cosmological observations.This paper aims at providing such models for stationary fields, which can be estimated from one observed sample only.
At thermal equilibrium, physical systems are usually characterized by the Gibbs probability distribution, also called Boltzmann distribution, that depends on the energy of the systems (Landau & Lifshitz 2013).For non-equilibrium systems, at a fixed time one may still specify the probability distribution of the field with a Gibbs energy, which is an effective Hamiltonian providing a compact representation of its statistics.Gibbs energy models can be defined as maximum entropy models conditioned by appropriate moments (Jaynes 1957).The main difficulty is to define and estimate the moments which specify these Gibbs energies.
For stationary fields, whose probability distributions are invariant to translation, moments are usually computed with a Fourier transform, which diagonalizes the covariance matrix of the field.The resulting covariance eigenvalues are the Fourier power spectrum.However, capturing non-Gaussian properties requires to go beyond second-order moments of the field.Third and fourth-order Fourier moments are called bispectrum and trispectrum.For a cubic -dimensional stationary field of length , the number of coefficients in the raw power spectrum, bispectrum and trispectrum are  (  ),  ( 2 ) and  ( 3 ) respectively.High-order moment estimators have high variance and are not robust, especially for non-Gaussian fields, because of potentially rare outliers which are amplified.It is thus very difficult to accurately estimate these high-order Fourier spectra from a few samples.Accurate estimations require to considerably reducing the number of moments and eliminating the amplification effect of high-order moments.
Local conservation laws for mass, energy, momentum, charge, etc. result in continuity equations or transport equations.The resulting probability distributions of the underlying processes thus are typically regular to deformations that approximate the local transport.These properties have motivated many researchers to use of a wavelet transform as opposed to a Fourier transform, which provides localized descriptors.Most statistical studies have concentrated on second-order and marginal wavelet moments (e.g., Bougeret et al. 1995;Vielva et al. 2004;Podesta 2009) which fail to capture important non-Gaussian properties of a field.Other studies (Ha et al. 2021) use wavelet operator for interpretation with application to cosmological parameter inference, but rely on a trained neural network model.
In recent years, new representations have been constructed by applying point-wise non-linear operators on the wavelet transforms of non-Gaussian fields to recover their high-order statistics.The scat- Steps to build a feasible model for a random field  from only one or a few realizations.We first build a low-dimension representation Φ( ) of the random field, which specifies a maximum entropy model.The representation Φ( ) is obtained by conducting the wavelet transform   and its modulus |  |, and then computing the means and covariance of all wavelet channels (  , |  |).Such a covariance matrix is further binned and sampled using wavelets to reduce its dimensionality, which is called the scattering spectra S ( ).Finally, These scattering spectra are renormalized and reduced in dimension by thresholding its Fourier coefficients along rotation and scale parameters Φ( ) =  S, making use of the regularity properties of the field.For many physical fields, this representation can be as small as only around ∼ 10 2 coefficients for a 256×256 field.
tering transform, for instance, is a representation that is built by cascading wavelet transforms and non-linear modulus (Mallat 2012;Bruna & Mallat 2013).This representation has been used in astrophysics and cosmology (Cheng & Ménard 2021), to study the interstellar medium (Allys et al. 2019;Saydjari et al. 2021), weaklensing fields (Cheng et al. 2020), galaxy surveys (Valogiannis & Dvorkin 2022), or radio observations (Greig et al. 2022).Other representations, which are built from covariances of phase harmonics of wavelet transforms (Mallat et al. 2020;Zhang & Mallat 2021), have also been used to model different astrophysical processes (Allys et al. 2020;Jeffrey et al. 2022;Régaldo-Saint Blancard et al. 2023).Such models, which can be built from a single image, have in turn enabled the development of new component separation methods (Regaldo-Saint Blancard et al. 2021;Delouis et al. 2022), which can be directly applied to observational data without any particular prior model of the components of a mixture (Auclair et al. 2023;Siahkoohi et al. 2023).
These models however suffer from a number of limitations: they are not very good at reproducing vortices or long thin filaments, and they require an important number of coefficients to capture dependencies between distant scales, as well as angular dependencies.Building on those previous works, reduced scattering covariance representations have been introduced, but only for time series, by leveraging scale invariance (Morel et al. 2022).In this paper, we present the scattering spectra, a low-dimensional representation that is able to efficiently describe a wide range of non-Gaussian processes encountered in physics.In particular, we show how it is possible to take into account the intrinsic regularity of physical fields to dramatically reduce the dimension of such representations.The first part of the paper presents maximum entropy models and scattering spectra statistics, as well as their dimensional reduction.The second part of the paper presents a quantitative validation of these models on various two-dimensional multiscale physical fields and discuss their limitations.
Notations:  * is the complex conjugate of a scalar .Ave  averages values indexed by  in a finite set.x [] is the Fourier transform of  [], whether  is a continuous variable in R  or belongs to finite periodic lattice.E{Φ()} is the expectation of Φ() according to the probability distribution () of a vector .log stands for base 2 logarithm.

Gibbs Energy of Stationary Fields
We review the properties of Gibbs energies resulting from maximum entropy models conditioned by moment values (Geman & Geman 1984;Zhu et al. 1997Zhu et al. , 1998)).We write  [] a field where the site index  belongs to a cubic -dimensional lattice of size .It results that  ∈ R   .
Assume that  ∈ R   has a probability density () and consider Gibbs energy models linearly parameterized by a vector  = {  } ≤ over a potential vector They define exponential probability models ) . (1) The model class is thus defined by the potential vector Φ(), which needs to be chosen appropriately.
If it exists, the maximum entropy distribution conditioned by E{Φ()} is a   0 which belongs to this model class.It has a maximum entropy In statistical physics,   0 is a macrocanonical model defined by a vector E{Φ()} of observables.One can verify that  0 also minimizes the Kullback-Liebler divergence within the class The main topic of the paper is to specify Φ() in order to define accurate maximum entropy models for large classes of physical fields, which can be estimated from a small number  of samples x .In this section, we suppose that  = 1.Reducing the model error given by (3) amounts to defining Φ which reduces the excess entropy of the model.This can be done by enriching Φ() and building very high-dimensional models.However, we must also take into account the empirical estimation error of E{Φ()} by Φ( x1 ), measured by E{ Φ() − E{Φ()}2 }.
In this paper, macrocanonical models are approximated by microcanonical models, which have a maximum entropy over a microcanonical set of width  > 0 Appendix A reviews a sampling algorithm for such model.It also explains how to extend the definition of Ω  for  > 1 samples x by replacing Φ( x1 ) by Ave  Φ( x ).If Φ() concentrates around E{Φ()} then the microcanonical model converges to the macrocanonical model when the system length  goes to ∞ and  goes to 0. The concentration of Φ() generally imposes that its dimension  is small relatively to the dimension   of .The choice of Φ() must thus incorporate a trade-off between the model error (3) and the distance between micro and macrocanonical distributions.

Fourier Polyspectra Potentials
Gaussian random fields are maximum entropy models conditioned on first and second-order moments.The potential vector Φ() is then an empirical estimator of first and second-order moments of .
For stationary fields, there is only one first-order moment E{ []} which can be estimated with an empirical average1 over : Similarly, the covariance matrix E{ [] [ ]} only depends on − , so only the diagonal coefficients in Fourier space are informative, which are called the power spectrum, The off-diagonal elements vanish because of phase cancellation under all possible translations, which means the second-order moments treat Fourier coefficients independently, and cannot describe relations or dependence between them.The diagonal elements, which can also be written as | x []| 2 , can be estimated from a single sample  by averaging | x [] | 2 over frequency bins that are large enough to reduce the estimator variance.A uniform binning and sampling along frequencies results in power spectrum estimators with  (  ) elements, so the Gaussian model is compact and feasible.However, the Gaussian random field model has limited power to describe complex structures.The majority of fields encountered in scientific research are not Gaussian.Non-Gaussianity usually means dependence between Fourier coefficients at different frequencies.The traditional way goes to higher orders moments of x, the polyspectra (Brillinger 1965), where phase cancellation implies that for stationary fields, only the following moments are informative, while other moments are zero.These polyspectra at order  > 2 capture dependence between  − 1 independent frequencies.As the leading term, the Fourier bispectrum specifies the non-zero thirdorder moments and has  ( 2 ) coefficients.However, bispectrum is usually not sufficient to characterize non-Gaussian fields.For example, it vanishes if the field distribution is symmetric () = (−).
One must then estimate fourth-order Fourier moments, the trispectrum, which has  ( 3 ) coefficients.There are two main problems for the polyspectra coefficients to become proper potential functions Φ() in the maximum entropy models.First, the number of coefficients increases sharply with the order.Second, high-order moments are not robust and difficult to estimate from a few realizations (Huber 1981).For random fields with a heavy tail distribution, which is ubiquitous in complex systems (Bak et al. 1987;Bouchaud & Georges 1990;Coles & Jones 1991;Kello et al. 2010;Sor 2017), higher order moments may not even exist.Those two problems are common for high-order moments and have been demonstrated in real-world applications (Dudok de Wit 2004;Lombardo et al. 2014).In the following two sections, we introduce modifications to this approach to solve those problems.

Wavelet Polyspectra
Many physical fields exhibit multiscale structures induced by nonlinear dynamics, which implies regularity of () in frequency.The wavelet transform groups Fourier frequencies by wide logarithmic bands, providing a natural way to compress the Fourier polyspectra.The compression not only reduces the model size but also improves estimator convergence.We use the wavelet transform to compute a compressed power spectrum estimate, as well as a reduced set of  (log 2 ) third and  (log 3 ) fourth order wavelet moments, allowing for efficient estimation of the polyspectra.

Wavelet Transform
A wavelet is a localized wave-form  [] for  ∈ R  which has a zero average ∫ R   []  = 0. We shall define complex-valued wavelets  [] = []    .where [] is a real window whose Fourier transform ĝ[] is centered at  = 0 so that ψ [] = ĝ[ − ] is localized in the neighborhood of the frequency .Fig. A1 shows  and ψ for a  = 2 dimensional Morlet wavelet described in appendix B. The wavelet transform is defined by rotating  [] with a rotation  in R  and by dilating it with dyadic scales 2  > 1.It defines Its Fourier transform is ψ which is centered at the frequency  and concentrated in a ball whose radius is proportional to 2 −  .
To decompose a field  [] defined over a grid of width , the wavelet is sampled on this grid.Wavelet coefficients are calculated as convolutions with periodic boundary conditions It measures the variations of  in a spatial neighborhood of  of length proportional to 2  , and it depends upon the values of x in a frequency neighborhood of  =  of length proportional to 2 −  .The scale 2  is limited to 1 ≤  ≤ , and for practical application to fields with a finite size , the choice of  is limited by  < log .Left part of Fig. 1 illustrates the wavelet transform of an image.The rotation  is chosen within a rotation group of cardinal , where  does not depend on .Wavelet coefficients need to be calculated for /2 rotations because  [, −] =  [, ] * for real fields.In  = 2 dimensions, the  rotations have an angle 2ℓ/, and we set  = 8 in all our numerical applications, which boils down to 4 different wavelet orientations.The total number of wavelet frequencies  is  =  (log ) 2 as opposed to   Fourier frequencies.
A wavelet transform is also stable and invertible if  satisfies a Littlewood-Paley condition, which requires an additional convolution with a low-pass scaling function  0 centered at the frequency  = 0.The specifications are detailed in appendix B.

Wavelet Power Spectrum
Given scaling regularity, one can compress the  (  ) power spectrum coefficients into  =  (log ) coefficients using a logarithmic binning defined by wavelets.This is obtained by averaging the power spectrum with weight functions as the Fourier transform of wavelets, which are band-pass windows, The limited number of wavelet power spectrum coefficients has reduced estimation variance.In fact, they are also the diagonal elements of the wavelet covariance matrix,  [, ] [, ] * = | [, ]| 2 , therefore an empirical estimation can also be written as an average over : Similar to the power spectrum, phase cancellation due to translation invariance means that the off-diagonal blocks i.e. the crosscorrelations between different wavelet frequency bands are nearly zero because the support of two wavelets ψ and ψ are almost disjoint, as illustrated in Fig. 2(a).

Selected 3rd and 4th Order Wavelet Moments
One may expect to compress the polyspectra in a similar manner with a wavelet transform, taking advantage of the regularities of the field probability distribution.However, it is non-trivial to logarithmically bin the polyspectra because more than one independent frequency is involved and the phase cancellation condition needs to be considered.
To solve this problem, let us revisit the phase cancellation of two frequency bands, which causes their correlation to be zero, for  ≠  .To create a non-zero correlation, we must realign the support of  [, ] and  [ ,  ] in Fourier space through nonlinear transforms.As shown in Fig. 2(b), we may apply a square modulus to one band (shown in blue) in the spatial domain, which recenters its frequency support at origin.Indeed, | ★   | 2 = ( ★   )( ★   ) * has a Fourier support twice as wide as that of  ★   , and will overlap with another wavelet band with lower frequency than .The transformed fields | ★   | 2 can be interpreted as maps of locally measured power spectra.Correlating this map with another wavelet band  ★   gives some third-order moments that are a priori non-zero.Furthermore, for wide classes of multiscale processes having regular power spectrum, it suffices to only keep the coefficients at  =  because of random phase fluctuation (see appendix B).For stationary random fields, they can be estimated with an empirical average over , to the ratio between the surface area of a -1-sphere and the volume of -1-ball, proportionally to Γ(/2 + 1/2)/Γ(/2).It results in an approximate scaling of  =  ( log ) when  is small and  ( Now we obtain a set of statistics characterizing the dependence of Fourier coefficients in two wavelet bands in a collective way, which are selected third-order moments.They can be interpreted as a logarithmic frequency binning of certain bispectrum coefficients.There are about  2  2 =  (log 2 ) such coefficients, which is a substantial compression compared to the  ( 2 ) full bispectrum coefficients.
Similarly, we consider the cross correlation between two wavelet bands both transformed by the square modulus operation and obtain a wavelet binning of fourth-order moments, For stationary fields, this covariance only depends on  −  .A further reduction of such a large covariance function is possible because its Fourier transform over  −  has two properties.First, it typically does not have higher frequency components than the initial wavelet transforms involved (see Fig. 2) as the phase fluctuations have been eliminated by the square modulus, and second, for fields with multiscale structures, it is regular and can be approximated with another logarithmic frequency binning.Thus, we can compress the large covariance function with a second wavelet transform, and estimate it by an empirical average over : where and the central frequencies of the second wavelets verifies || ≥ | | > ||.There are about  3  3 =  (log 3 ) such coefficients, which is also a substantial compression compared to the  ( 3 ) full trispectrum coefficients.

Scattering Spectra
In general, the estimation of high-order moments has a high variance because high-order polynomials amplify the effect of outliers.A scattering approach (Mallat 2012;Bruna & Mallat 2013;Cheng & Ménard 2021) reduces the variance of these estimators by replacing || 2 by ||.The resulting spectra only depend on the mean and covariance matrix of (, ||), which are low-order transforms of the original field .
Local statistics of wavelet modulus have been studied to analyze properties of image textures (Portilla & Simoncelli 2000).Their mathematical properties have been analyzed to capture non-Gaussian characteristics of random fields (Mallat et al. 2020;Zhang & Mallat 2021) in relation to scattering moments (Mallat 2012;Bruna & Mallat 2013).Scattering spectra have been defined on one-dimensional time-series (Morel et al. 2022), from the joint covariance of a wavelet transform and its modulus: ( , ||).We extend it to fields of arbitrary dimension  and length , in relation to Fourier high-order moments, and define models of dimension  (log3 ).

First and second wavelet moments, sparsity
For non-Gaussian fields , wavelet coefficients  [, ] define fields which are often sparse (Olshausen & Field 1996;Stephane 1999).This is a non-Gaussian property that can be captured by first-order wavelet moments E{| [, ]|}.If  is a Gaussian random field then  [, ] remains Gaussian but complex-valued so, and we have 4 .This ratio decreases when the sparsity of  [, ] increases.The expected value of || is estimated by and the ratio is calculated with the second-order wavelet spectrum estimator

Cross-Spectra between Scattering Channels
A scattering transform is computed by cascading modulus of wavelet coefficients and wavelet transforms (Mallat 2012;Bruna & Mallat 2013).Let us replace || 2 by || in the selected third and fourth-order wavelet moments described in the previous section.
Replacing || 2 by || in the fourth order wavelet moments (11) amounts to estimating the covariance matrix of wavelet modulus fields ||.As the  −  dependency of this covariance can also be characterized by a second wavelet transform, this amounts in turn to estimate the covariance of scattering transforms for || ≥ | | ≥ ||.It provides a wavelet spectral estimation of the covariance of ||.
Combining the moment estimators of Eqs.(12,13,14,15) defines a vector of scattering spectra It provides a mean and covariance estimation of the joint wavelet and wavelet modulus vectors (, ||).It resembles the second, third, and fourth-order Fourier spectra but has much fewer coefficients and better information concentration.Considering the conditions satisfied by ,  , and , the exact dimension of () is  +  2  ( − 1)/8 +  3  ( 2 − 1)/48, of the order  ( 3 ).

Renormalization
Scattering spectra coefficients must often be renormalized to improve the sampling of maximum entropy models.Indeed, multiscale random processes often have a power spectrum that has a power law decay E{| x []| 2 } ∼ | | − over a wide range of frequencies, long-range correlations corresponding to a strong decay from large to small scales.The wavelet spectrum also has a power-law decay This means that if we build a maximum entropy model with Φ() = () then the coordinate of Φ() of low-frequencies  have a much larger amplitude and variance than at high frequencies.The microcanonical model is then dominated by low frequencies and is unable to constrain high-frequency moments.
The same issue appears when computing the  0 parameters of a macrocanonical model defined in ( 2), for which it has been shown that renormalizing to 1 the variance of wavelet coefficients at all scales avoid numerical instabilities (Marchand et al. 2022) 3 .We renormalize the scattering spectra by the variance of wavelet coefficients,  2 [] = Ave   2 ( x ) [], which can be estimated from a few samples.The renormalized Scattering Spectra are .
The microcanonical models proposed in this paper are built from these renormalized statistics and/or their reduced version described below.

Dimensionality reduction for physical fields
Though much smaller than the polyspectra representation, the scattering spectra S representation still has a large size.Assuming isotropy and scale invariance of the field , a first-dimensional reduction can be performed that relies on the equivariance properties of scattering spectra with respect to rotation and scaling (see appendix C).However, such invariances cannot be assumed in general.In this section, we propose to construct a low-dimensional representation by only assuming regularity under rotation or scaling of the scales involved in the scattering spectra representation.A simplified version of such a dimensional reduction has been introduced in (Allys et al. 2019).
We refer the reader to appendix D for technical details.
The goal of the reduction is to approximate the covariance coefficients S3 and S4 , the most numerous, using only a few coefficients.This can be seen as a covariance matrix estimation problem.To do so, we first use a linear transform to sparsify the covariance matrix and then perform a threshold clipping on the coefficients to reduce the representation.We consider a linear transform  S = ( S1 , S2 ,  S3 ,  S4 ) with a pre-determined linear transform  which stands for a 2D or 3D Fourier transform along all orientations, as well as a 1D cosine transform along scales, for S3 and S4 .For fields with statistical isotropy or self-similarity, all harmonics related to the action of global rotation and scaling on the field  should be consistent with zero, except for the zeroth harmonic.For general physical fields, we expect the statistics S() to have regular variations to the action of rotation or scaling of the different scales involved in its computation, which implies that its Fourier harmonics  S() have a fast decay away from the 0-th harmonic and  S() is a sparse representation.
Thresholding on a sparse representation is widely used in image processing for compression (Chang et al. 2000).We use threshold clipping on the sparse representation  S to significantly reduce the size of the scattering spectra.Furthermore, when empirically estimating large but sparse covariance matrices such as , thresholding provides Stein estimators (Stein 1956) which have lower variance and are consistent(e.g., Donoho & Johnstone 1994;Bickel & Levina 2008;Cai & Liu 2011;Fan et al. 2013).As S1 or S2 are already small, we keep all of their coefficients.
There are different strategies available to set the threshold for clipping.We adopt a simple strategy which keeps those coefficients with ( S) > 2( S), where ( S) and ( S) are the means and standard deviations of individual coefficients of  S.These adaptive thresholding estimators achieve a higher rate of convergence and are easy to implement (Cai & Liu 2011).With multiple realizations from simulations, ( S) and ( S) can be estimated directly.In the case where only a single sample field is available, ( S) can be estimated from different patches of that sample field (e.g., Sherman 2018).We call  S the coefficients after thresholding projection: The compact yet informative set of scattering spectra  S is the representation Φ() =  S() proposed in this paper to construct maximum entropy models.

NUMERICAL RESULTS
We have introduced maximum entropy models based on small subsets of  (log 3 ) scattering spectra moments S and projected moments  S, claiming that it can provide accurate models of large classes of multiscale physical fields, and reproduce  ( 3 ) power spectrum, bispectrum and trispectrum Fourier moments.This section provides a numerical justification of this claim with five types of 2D physical fields from realistic simulations.In order to reduce the variance of the validation statistics, we consider in this section a model estimated on several realizations of a field.However, our model also produces convincing realizations when estimated on a single realization (see Fig. B1 for a visual assessment).

Dataset of Physical Fields
We use five two-dimensional physical fields to test the maximum entropy models.The five fields are chosen to cover a range of properties in terms of scale dependence, anisotropy, sparsity, and morphology: (D) Magnetic turbulence: column density of 3D isothermal magnetic-hydrodynamic (MHD) turbulent simulations (Allys et al. 2019).The field is anisotropic due to a mean magnetic field in the horizontal direction.(E) Anisotropic turbulence: two-dimensional slices of a set of 3D turbulence simulations (Li et al. 2008;Perlman et al. 2007).To create anisotropy, we have squeezed the fields along the vertical direction.
These simulations are sampled on a grid of 256×256 pixels with periodic boundary conditions4 and normalized to have zero mean and unity standard deviation, respectively.Samples of each field are displayed in the first row of Fig. 3. To clearly show the morphology of small-scale structures, we zoom in to a 128×128 region.

Model description and visual validation
We fit our maximum entropy model using wavelet polyspectra and scattering spectra, respectively, with the following constraint, where the second average is computed on an ensemble of 100 realizations x for each physical simulation (for field D we use only 20 realizations due to the availability of simulations), and the field generation is performed simultaneously for 10 fields   , making our microcanonical model closer to its macrocanonical limit.The microcanonical sampling algorithm is described in appendix A.
Examples of field generation results are given in Fig. 3.The second row shows samples generated based on the high-order normalized wavelet moments are defined similarly to S in (17).For the choice of wavelets, we use J=7 dyadic scales, and we set  = 8 which samples 4 orientations within , resulting in dim M = 11 677 coefficients for M. The third row in Fig. 3 shows results from a reduced set Φ() =  M (), which is a 2 Fourier thresholded representation of M defined in exactly the same way as  S in (18).The thresholding yields dim  M = 147, 286, 547, 1708, 926 for fields A-E, respectively.A visual check shows that these models fail to recover all morphological properties in our examples especially when a thresholding reduction is applied.This issue is a manifestation of the numerical instability of high-order moments.
In the fourth row, we present sample fields modeled with the scattering spectra  with dim  S = 11 705 for J=7 and R=8.A visual check reveals its ability to restore coherent spatial structures including clumps, filaments, curvy structures, etc.The low-order nature and numerical stability of  also significantly fasten the sampling compared to the high-order moments M (200 vs. 800 steps to converge).The last row shows sample fields modeled by a much smaller set , which has dim  S = 204, 364, 489, 615, 304 coefficients for fields A-E, respectively.This model is ∼ 10 2 times smaller, while generating samples visually indistinguishable from the full set model with Φ() = ().In addition, the ratio between the dimensionality of the field dim  =   (the number of pixels) and the model dim Φ is more than 100.

Statistical Validation
We now quantify the consistency between the scattering spectra models and the original fields using a set of validation statistics  () defined below, including marginal PDF, structure functions   , power spectrum , and normalized bispectrum B and trispectrum T. The validation statistics are shown in Figs. 3 and 4, where black curves represent the expected value  original of these statistics, estimated from 100 realizations x of the original simulated fields (except for field D for which we have only 20 realizations).Gray regions around the black curves represent the standard deviations  original of those statistics estimated on the original fields.Blue curves are statistics  S,model estimated on fields modeled with .Similarly,   S model are estimated on fields modeled with the reduced set .Both these averages are estimated from the 10 fields simultaneously sampled from the corresponding microcanonical models.

Validation statistics
The marginal probability distribution function (PDF) is measured as the histogram of sample fields and shown in Fig. 3.It averages out all spatial information and keeps only the overall asymmetry and sparsity properties of the field.The marginal information is not explicitly encoded in the scattering spectra, but for all the five physical fields we examine here, it is recovered even with the reduced model  S, where only ∼ 10 2 scattering spectra coefficients are used.
Given that the high dimensionality of the full set of polyspectra coefficients, as well as the computational cost of estimating them properly, we adopt an isotropic shell binning for the power spectrum, bispectrum, and trispectrum.Although this reduces the number of coefficients as well as their variance, working with isotropic statistics prevents the characterization of anisotropic features, for instance in fields D and E, unlike with scattering spectra.Validation results with these isotropic polyspectra are given in Fig. 4.
The shell binning is defined as follow.We first divide the Fourier space into 10 annuli with the frequencies linearly spaced from 0 to 0.4 cycles/pixel.Then, we average the power and poly spectra coefficients coming from the same annulus combinations.For instance, the power spectrum yields: To decorrelate the information from the power spectrum and higher orders, we normalized the binned bi-and tri-spectra by []: , where the   -dimensional wave-vectors are respectively averaged in the  th  frequency annuli, and satisfy    = 0. To clearly reveal the diversity of different type of physical fields, the trispectrum T coefficients shown in Fig. 4 are subtracted by the reference value of Gaussian white noise, evaluated numerically on 1000 independent realizations.Details about the numbers and the ordering of B and T are given in appendix E.
In Fig. 4 we also show the validation with structure functions, which are -th order moments of the field increments as a function of the lag Initially proposed by Kolmogorov for the study of turbulent flows (Kolmogorov 1941), they are widely used to analyze non-Gaussian properties of multiscale processes (Jaffard 2004).

Comparison between original and modeled fields.
We quantify the discrepancy between the model and original field distributions by the outlier fraction of validation statistics outside the 2 range, For each of the five types of fields, we observe the following fractions.The binned power spectrum has fractions of : 0%, 0%, 20%, 0%, 0% for the models using all S statistics and 0%, 10%, 40%, 10%, 0% for the thresholding models with  S. The power spectrum deviation of field C is likely caused by the longer convergence steps required by smooth fields, as our generative models start from white noise with strong small-scale fluctuations.Indeed increasing the steps to 800 reduces the outlier fraction of the  S model to 10%.For B and T, the outlier fractions are all below 5% except for the models of field A, where the bispectrum coefficients have 13% of outliers.Those outliers all have the smallest scale involved, and disappear if the high-frequency cut is moved from 0.4 to 0.35 cycles/pixel.The low fractions demonstrate consistency between our maximum entropy models and ensembles of the original physical fields.
For field A, a similar deviation is also observed in high-order structure functions.For this field, it can be seen from Fig. 4 that even though many coefficients are not defined as outliers, they all tend to have a lower value than the original ones.This effect may originate from the log-normal tail of the cosmic density field (Coles & Jones 1991), whose Gibbs potential includes terms in the form of log , in contrast to the form of || in scattering covariance or   in high-order statistics.However, regardless of this difficulty, these outliers are all still within a 3 range, demonstrating that the scattering spectra provide a good approximation though not exact model for fields with such heavy tails.
The marginal PDF, structure functions, power spectrum and polyspectra probe different aspects of the random field ().The polyspectra especially probe a huge variety of feature configurations.For all the validation statistics, we observe general agreement between the model and original fields.Such an agreement is a nontrivial success of the scattering spectra model, as those statistics are not generically constrained by the scattering spectra for arbitrary random fields.They indeed significantly differ from the scattering spectra in the way they combine spatial information at different frequencies and in the non-linear operation adopted.The agreement implies, as we have argued, that symmetry and regularity can be used as strong inductive bias for physical fields and the scattering spectra, with those priors build-in, can efficiently and robustly model physical fields.

Visual Interpretation of Scattering Spectra Coefficients
The key advantage of the scattering spectra compared to usual convolutional neural networks is their structured nature: their computation corresponds to the combination of known scales and orientations in a fixed way.Beyond the limited number of symmetries, the structured nature of the scattering spectra allows us to both quantify and interpret the morphology of structures, which is one of the original goals to design these statistics.
The values of scattering spectra can be shown directly (see Fig. C1) to analyze non-Gaussian properties of the field.Moreover, the meaning of its coefficients can also be visualized through our maximum entropy generative models.As one gradually changes the value of some summary statistics, the morphology of structures in the generated fields also changes.A similar exploration for a smaller set of scattering transform coefficients has been explored in Cheng & Ménard (2021), and we show such results with the much more expressive scattering spectra coefficients in Fig 5 .Such exploration using synthesis is also similar to the feature visualization efforts for convolutional neural networks (Olah et al. 2017).
The central panel is a realization of field B from physical simulations.The other four panels are generated fields with two collective modifications of the scattering spectra: the vertical direction shows the effect of multiplying all S3 and S4 coefficients by a factor of 1/3 or 3.It indicates that the amplitude of S3 and S4 controls the overall non-Gaussian properties of the field and in particular the sparsity of its structures.The horizontal direction corresponds to adjusting the orientation dependence.We set the coefficients with parallel wavelet configurations (i.e., S3 ) as references and keep them unchanged.Then, we make the difference from other coefficients to those references to be 2 times or -2 times the original difference.Visually, it controls whether structures are more point-like or more curvy-like in the field.In this experiment, the generated field is initialized with the original field instead of white noise, in order to clearly show the correspondence between the field structure and scattering spectra coefficients.

Application to Identifying Symmetry
As an expressive representation whose coefficients are equivariant under standard group transformation, the scattering spectra can also be used to detect and identify the various statistical invariances commonly present in physical fields.Besides the aforementioned rotation and scaling invariance, more can also be included, such as the flipping of coordinate or field values.
The simplest way to check asymmetry to a transformation like rotation or flip is to check if the scattering spectra  are changed after applying such a transform.A more sophisticated way that can also quantify partial symmetries is to linearly decompose S into symmetric and asymmetric parts and then compute the fraction of asymmetric coefficients surviving the thresholding reduction.We further normalize this fraction by that in the full set: .
When it is zero, the random field () should be invariant to the transform up to the expressivity of our representation.For the five random fields analyzed in this study, we measure their asymmetry indices with respect to rotation and scaling.The corresponding anisotropy and scale dependence indices are (A) 0, 0.16 ; (B) 0, 0.53; (C) 0, 0.66; (D) 0.32, 0.45; (E) 0.28, 0.29.As expected, the cosmic lensing field (field A), which consists of haloes at all scales and strengths, is closest to isotropic and scale-free.The cosmic web (B) and 2D turbulence (C) fields are isotropic but have particular physical scales above which the field becomes Gaussian, so they are not scale-free.The last two turbulence fields have anisotropic physical input, but the latter largely probes the 'inertial regime' of turbulence, which is scale-free.

Limitations
While a broad range of physical fields satisfy the implicit priors of the scattering covariance, one does expect regimes for which the description will not be appropriate.The so-called  4 field in physics comes as a first problematic example.It is the maximum entropy field under the power spectrum and pointwise fourth-order moment  4 constraints, but this characterization is unstable to specify a non-convex pdf which is a pointwise property as opposed to the delocalized Fourier moments and it is highly unstable at critical points (Marchand et al. 2022).The first column in Fig. 6 shows an original  4 field at its critical temperature and that generated from the full set of scattering covariance.In contrast to previous examples, this type of field is not successfully reproduced.
On the other hand, when built based on one example field  1 and generating only one realization x1 (i.e., in Eq. 19 both  and  are 1), our model has a risk of over-fitting: it almost exactly copies the original field with an arbitrary translation and does not provide enough randomness.It can also be seen as a transition from generative modeling regime into a coding regime.This is related to the fact that for maximum entropy models, when the number of constraints amounts to a considerable fraction of the number of total degree of freedom, the microcanonical distribution deviates significantly from the macrocanonical distribution, and has a much lower entropy.The middle panel of Fig. 6 illustrate this effect, where the relative position of triangles of the modeled field is exactly copied from the original field.It happens only when the field is sparse, and when the full set S is used.This problem can be avoided by increasing the number of input fields or generated fields, or an early stop in the microcanonical sampling.
For physical fields with multi-scale structures, it is expected that the distribution function () does not change much under a slight deformation.When modeling such fields, it is important to have a representation that has the same property.Being built from wavelet decomposition and contracting operator, the scattering spectra also linearize small deformation in the field space, which plays an important role in lowering its variance (see (Bruna & Mallat 2013)).However, when modeling structured fields whose distribution functions are not regular under deformation, this means that the generative model will simply produce structures that are "close enough" up to small deformations.This typical type of failure is shown in the third example of Fig. 6.

CONCLUSION
We build maximum entropy models for non-Gaussian random fields based on the scattering spectra statistics.Our models provide a lowdimensional structured representation that captures key properties encountered in a wide range of stationary physical fields, namely: (i) stability to deformations as a result of local conservation laws in Physics for mass, energy, momentum, charge, etc; (ii) invariance and regularity to rotation and scaling; (iii) scale interactions typically not described by high-order statistics; Those are the priors included in the scattering spectra.
Our models provide a practical tool for generating mock fields based on some example physical fields.In sharp contrast to neural network models, our representation has the key advantage of being interpretable and can be estimated on a few realizations.This is crucial in Physics where generating fields in experiments or simulations is costly or when non-stationarity limits the amount of clean recorded data.Our proposed approach enables a new range of data/simulation analyses (e.g.Regaldo-Saint Blancard et al. 2021;Delouis et al. 2022), involving extensions to the modeling of cross-regularities when multiple channels are available (e.g.Régaldo-Saint Blancard et al. 2023).

APPENDIX A: MICROCANONICAL SAMPLING
Given  observed samples  1 , . . .,   of a field, with possibly  = 1, the microcanonical ensemble given (4) can be extended as follow: Microcanonical models are maximum entropy distributions over Ω  , which have a uniform distribution over this ensemble.Increasing the number of samples  reduces the variance of Ave  Φ(  ) which concentrates around E{Φ()}.This reduces the information about a specific realization which is contained in Ave  Φ(  ), thus limiting over-fitting.
Sampling from the microcanonical model amounts to drawing a realization from a uniform distribution in Ω  .We approximate this sampling with a gradient descent algorithm studied in (Bruna & Mallat 2019).This algorithm progressively transports a white Gaussian noise distribution, which has a higher entropy than the microcanonical model, into distributions supported in Ω  .This is done with a gradient descent on ℓ( 1 , . . .,   ) = Ave  Φ(  ) − Ave  Φ( x ) 2 , where the   are initialized as independent realizations of white noises.At each iteration, the   are updated with the L-BFGS-B algorithm, which is a quasi-Newton method that uses an estimate of the Hessian matrix.In practice, we perform 200 gradient descent steps which yield a typical error  ≈ 10 −4 .
It is proved in (Bruna & Mallat 2019) that this algorithm converges to a distribution that has the same symmetries as Φ(), similarly to the microcanonical one.However, it has been shown that this algorithm recovers a maximum entropy distribution in Ω  only under appropriate conditions and that such gradient descent models may differ, in general, from maximum entropy ones.Nevertheless, these algorithms provide powerful sampling methods to approximate large classes of high-dimensional stationary processes, while being much faster and computationally tractable than alternative MCMC algorithms.
depends on  −  and have a fast decay when the power spectrum () is regular.Thus, even if dependencies across separate scales may exist, they are not captured by correlation.
Taking the modulus of wavelet coefficients removes complex phase oscillations and thus recenter the frequency support of  [, ].Indeed, the power spectrum   () of  ★   is mostly supported in a ball  −  ≤ 2 −   −1 which does not overlap with the Fourier support of the power spectrum   () of ★  .Taking a modulus on  ★   eliminates the phase which oscillates at the central frequency  .As a consequence, the power spectrum of | ★   | is centered at  = 0 and its energy is mostly concentrated in  ≤ 2 −   −1 which now may overlap with the support of   () [𝜆] as can be seen in Fig. 2. The power spectra of | [, ]| and | [,  ] |, both centered at zero, also overlap.
We now justify taking  =  in order 3 moments given by (10).The cross spectrum  , () between  [, ] and | [,  ] | is assumed regular for the fields considered in this paper.In that case one can approximate such cross-spectrum using wavelets, which gives the moments E{ [, , ]  ||[,  , ]}.However, the lefthand-side  [, , ] is negligible when  ≠  because Fourier support of wavelets   and   barely overlap.The resulting coefficients For the physical fields studied in this paper, such coefficients are non-zero, thus revealing their non-Gaussianity Fig. C1.

APPENDIX C: EQUIVARIANCE AND INVARIANCE TO ROTATIONS AND SCALING
The scattering spectra are computed from wavelet transforms, which are equivariant to rotations and scalings.We show that scattering spectra inherit these equivariance properties.If () is isotropic or self-similar, then one can build isotropic or self-similar maximum entropy models by averaging renormalized scattering spectra over rotations or scales, which reduces both the variance and dimensionality of S.
To avoid discretization and boundary issues for rotations and scaling, we consider fields  [] defined over continuous variables  ∈ R  , and establish the mathematical results in this framework.For this purpose, the sum in the wavelet transform defined in ( 8) is replaced by an integral over R  .Wavelets are dilated by 2  for  ∈ Z and rotated by  in a rotation group  of cardinal .In dimension  = 2, these rotations have an angle 2ℓ/.Power-spectrum  2 and sparsity factors S1 are averaged along all angles (amount to taking the 0-th angle Fourier harmonic).We only show the 0th angle Fourier harmonic and 0-th scale Fourier harmonic for order 3 and order 4-moment estimators S3 and S4 .Thus, the quantities that are shown are invariant to the rotation of the field, and the last two rows ( S3 , S4 are furthermore invariant to scaling).Non-zero coefficients S3 show that the cosmic lensing and cosmic web fields are not invariant to sign flip.This is due to the presence of high positive peaks on the former and filaments on the latter.The large amplitude of envelope coefficients S4 on the last 2 fields indicate long-range spatial dependencies as evidenced by the presence of structures at the level of the map.transform calculation.However, the same approach applies to noncommutative groups  of rotations in R  for  > 2, with their Fourier transform.Each wavelet frequency is defined in (7) by  = 2 −   ℓ , where  ℓ is a rotation of angle 2ℓ/.To guaranty that the scattering spectra frequencies satisfy || ≤ | | < ||, we write If S() has regular variations as a function of rotations then its three-dimensional Fourier transform along the (ℓ 1 , ℓ 2 , ℓ 3 ) has coefficients of negligible amplitude at high frequencies, which can thus be eliminated.One can also take advantage of regularities along scales.Since 1 ≤  1 ≤  varies on an interval without periodicity, the Fourier transform is replaced by a cosine transform along  1 for  and  fixed.We could also perform a cosine transform along the scale shift  and , but this is not done in numerical applications because their range of variations is small and -dependent.The Fourier transforms along  1 is however sufficient to identify scale-invariance, since one then expects S to only depend on  and , see appendix C. We write  S() the Fourier transform of S() along (ℓ 1 , ℓ 2 , ℓ 3 ) and its cosine transform along  1 .
Since  is unitary, it preserves the estimator variance: Ideally, the estimation error of E{ S()} is reduced by eliminating its coefficients whose squared amplitude is smaller than the variance of the empirical estimation error.It amounts to suppressing all coefficients having a variance that is larger than the bias resulting from their elimination.However, we can not implement this optimal "oracle" decision because we do not know E{ S()}.In this paper, we instead apply an approximate thresholding algorithm, which eliminates small amplitude coefficients of S() below a threshold proportional to their standard deviation, as discussed in the main text.This thresholding algorithm is adaptive and the selected coefficients vary from one process to another.For each process studied, an ensemble of between 20 to 100 samples {  } were used to empirically estimate the average and variance of  S, called ( S) and ( S).If () is isotropic or self-similar then we expect that  is a lowfrequency projector along global rotations (which act similarly on all   coordinates) or scalings (which act on ), which corresponds to the averages described in (C7) and (C8).The Fourier projection  is however much more general and can adapt to unknown regularities of () along rotations and scales.
way.The frequency annuli are labeled by from small to large | |.
To remove redundant coefficients, we require  1 ≤  2 ≤  3 (≤  4 ) and order them first by  1 in increasing order; when two binning configurations have the same  1 , they are then ordered by  2 and so on.

Figure 1 .
Figure 1.Steps to build a feasible model for a random field  from only one or a few realizations.We first build a low-dimension representation Φ( ) of the random field, which specifies a maximum entropy model.The representation Φ( ) is obtained by conducting the wavelet transform   and its modulus |  |, and then computing the means and covariance of all wavelet channels (  , |  |).Such a covariance matrix is further binned and sampled using wavelets to reduce its dimensionality, which is called the scattering spectra S ( ).Finally, These scattering spectra are renormalized and reduced in dimension by thresholding its Fourier coefficients along rotation and scale parameters Φ( ) =  S, making use of the regularity properties of the field.For many physical fields, this representation can be as small as only around ∼ 10 2 coefficients for a 256×256 field.

( A )
Cosmic lensing : simulated convergence maps of gravitational lensing effects induced by the cosmological matter density fluctuations(Matilla et al. 2016;Gupta et al. 2018).(B) Dark matter: logarithm of 2D slices of the 3D large-scale distribution of dark matter in the Universe(Villaescusa-Navarro et al. 2020).(C) 2D turbulence: turbulence vorticity fields of incompressible 2D fluid stirred at the scale around 32 pixels, simulated from 2D Navier-Stokes equations(Schneider et al. 2006).

Figure 3 .
Figure3.Visual comparison of realistic fields and those sampled from maximum entropy models based on wavelet higher-order moments M and wavelet scattering spectra S statistics.The first row shows five example fields from physical simulations of cosmic lensing, cosmic web, 2D turbulence, magnetic turbulence, and squeezed turbulence.The second and third rows show syntheses based on the selected high-order wavelet statistics estimated from 100 realizations.They are obtained from a microcanonical sampling with 800 steps.The fourth and fifth rows show similar syntheses based on the scattering spectra statistics, with only 200 steps of the sampling run.This figure shows visually that the scattering spectra can model well the statistical properties of morphology in many physical fields, while the high-order statistics either fail to do so or converge at a much slower rate.To clearly show the morphology of structures at small scales, we show a zoom-in of 128 by 128 pixels regions.Finally, to quantitatively validate the goodness of the scattering model, we show the marginal PDF (histogram) comparison in the last row.1-15(2023)

Figure 4 .
Figure 4. Validation of the scattering maximum models for the five physical fields A-E by various test statistics.The curves for field E represent the original statistics and those for A-D are shifted upwards by an offset.In general, our scattering spectra models well reproduce the validation statistics of the five physical fields.

Figure 5 .
Figure 5. Visual interpretation of the scattering spectra.The central is one realization of field B in physical simulations.The other four panels are generated fields with two simple collective modifications of the scattering spectra coefficients.

Figure 6 .
Figure 6.Example of failures and applications typical physical fields.The modeled field of the central panel has been recentered for easier comparison with the original ones.

Figure B1 .
Figure B1.Visual assessment of our model based on S 11 641 coefficients estimated on a single realization (top).Generated fields (bottom) show very good visual quality.

Figure C1 .
Figure C1.Visualization of Scattering Spectra S different physical fields.Power-spectrum  2 and sparsity factors S1 are averaged along all angles (amount to taking the 0-th angle Fourier harmonic).We only show the 0th angle Fourier harmonic and 0-th scale Fourier harmonic for order 3 and order 4-moment estimators S3 and S4 .Thus, the quantities that are shown are invariant to the rotation of the field, and the last two rows ( S3 , S4 are furthermore invariant to scaling).Non-zero coefficients S3 show that the cosmic lensing and cosmic web fields are not invariant to sign flip.This is due to the presence of high positive peaks on the former and filaments on the latter.The large amplitude of envelope coefficients S4 on the last 2 fields indicate long-range spatial dependencies as evidenced by the presence of structures at the level of the map.
2 S = E{ Ave   S(  ) − E{ S()} 2 }. (D1) The coefficients which have been kept are those that individually verify ( S) > 2( S).A projected scattering spectraΦ() =  S()is computed with a linear Fourier projection  which eliminates all coefficients of  S() corresponding to coefficients of Ave   S(  ) below their threshold.The efficiency of this projected scattering is the variance reduction ratio  2  S / 2 S with  2  S = E{ Ave   S(  ) − E{ S()} 2 }. (D2)