Pierre Chainais
email: pchainai@isima.fr

Véronique Delouille
email: v.delouille@oma.be

Jean-François Hochedez
email: hochedez@oma.be

Modeling images of the Quiet Sun in the extreme ultra-violet

Keywords: solar physics, multifractal analysis, wavelets, multiscaling processes

We address the statistical modeling of solar images provided by the Extreme ultraviolet Imaging Telescope (EIT) onboard the Solar and Heliospheric Observatory (SoHO, a joint ESA/NASA mission). We focus in particular on the less structured regions, the "Quiet Sun". We first review on a brief historical viewpoint on multifractal processes for physical modeling. Then we present a multifractal analysis of Quiet Sun images. Our aim is to identify a model that would permit to simulate images that are similar to real ones, and to use the scale invariance property to obtain artificial images at any finer resolution. We compare various families of models including infinitely divisible cascades and fractional stable fields that permit to synthesize images that are statistically similar to Quiet Sun images. This modeling will assist in promoting forthcoming high resolution observations by analysing sub-pixel variability in today's solar corona images.

INTRODUCTION

The complex dynamics of the Sun obeys an approximately 11 years periodic cycle. During the most active phases, one observes very large magnetic structures with their specific geometry and dynamics. Such organized structures are called 'active regions', see figure 1(a). They are the subject of detailed studies. [START_REF] Moses | EIT observations of the extreme ultraviolet sun[END_REF][START_REF] Schrijver | A new view of the solar outer atmosphere by the transition region and coronal explorer[END_REF] During more quiet periods, such 'active regions' are less numerous ; in the bandwidth centered on λ = 19.5 nm, one observes large regions which are much less organized, in appearance at least. These so-called 'Quiet Sun' areas [START_REF] Berghmans | Quiet Sun EUV transient brightenings and turbulence. A panoramic view by EIT on board SOHO[END_REF][START_REF] Pauluhn | Statistics of quiet Sun extreme ultraviolet intensities[END_REF] look like some turbulent background, see figure 1(b), but their precise nature is still elusive. Despite its sounding peaceful, this name certainly hides a much more lively small scale reality. Small scales remain out of reach due to a lack of resolution of present instruments.

The Extreme ultraviolet Imaging Telescope (EIT) [START_REF] Delaboudinière | EIT: Extreme ultraviolet imaging telescope for the soho mission[END_REF] onboard the SoHO (Solar and Heliospheric Observatory) mission jointly launched by the European Space Agency (ESA) and the National AeroSpace Agency (NASA) provides images of the Sun in the extreme ultraviolet bandwidth, at the wavelength λ = 19.5 nm in particular. The resolution is about (1800km) 2 /pixel. This is rather rough compared to the smallest physical scales which are supposed to be as small as 100 m. Therefore the quality of observations is limited and we do not know really how far they are limited. We suggest to take benefit from the scale invariance property observed on these images to extrapolate the sub-pixel structure, at least in a statistical sense. Indeed, our aim is not to extract sub-pixel information as superresolution would. Our aim is to propose a statistical model of what sub-pixel information may be, given 'badly resolved' current observations. We expect that sub-pixel extrapolation based on scale invariance provide us with a model for Quiet Sun images at any resolution. This model should feature some statistical properties of images that are currently available. The idea is to consider these images as (disordered) stochastic textures. Even though we do not pretend that such a model will precisely take into account the whole physics of the system, it will not consist of a simple 'texture modeling'. Indeed, our approach relies on the use of stochastic multifractals models, the definition of which was initially motivated by the modeling of turbulent flows. [START_REF] Frisch | Turbulence. The legacy of A. Kolmogorov[END_REF][START_REF] Chainais | Multidimensional infinitely divisible cascades. application to the modelling of intermittency turbulence[END_REF] From a more general viewpoint, such an approach may be interesting when observing a disordered turbulent system with a resolution much bigger than the sallest physical scale (e.g., clouds in meteorology, star dusts in astrophysics). An essential motivation of this work is the anticipation of the quality of future highly resolved observations. [START_REF] Hochedez | Wide bandgap EUV and VUV imagers for the Solar Orbiter[END_REF] Physicists are studying the opportunity to reach resolutions of the order of (80 km) 2 /pixel, in place of the current (1800 km) 2 of EIT. The scientific interest of such expensive missions has to be argued. The CCD sensors used in Solar Physics are essentially photon-counting instruments. The arrival of photon can be modelled by a Poisson process. Among others, the problem of high resolution observations is the decrease in the signal to noise ratio as the resolution becomes finer. When the average number of incident photons becomes too small, the measurement is no more reliable. As a consequence, it is difficult to predict whether one may hope to gain some reliable detailed information by using a better resolution since we do not know a priori how the intensity is really distributed inside a presently 'badly resolved' pixel. The main thrust of this work is that we may propose an extrapolation (in the statistical sense) of the structure of images below the available pixel scale thanks to the use of scale invariant models. Eventually, we may be able to define a criterion to define regions where a high resolution might be relevant: regions with sufficiently high signal to noise ratio and probable presence of sub-pixel details.

The images synthesized by using our model will also permit to calibrate some tools such as Velociraptor 9 that was designed to extract the movements at the surface of the Sun thanks to an optical flow technique. For instance, we will be able to test the limits of such a tool when aiming at the extraction of movements much smaller than 1 pixel between to successive images * .

The paper is organized as follows. Section 2 recalls on the roots of multifractal analysis and multiplicative cascades for the modeling of physical systems. Section 3 recalls how wavelet are being used in the multifractal analysis of images. Section 4 presents our method to identify the parameters of stochastic multifractal models that capture a set of statistical properties of Quiet Sun images. Finally we discuss the validity and potential utility of our approach.

A BRIEF HISTORY OF MULTIFRACTAL MODELS FOR PHYSICS

From turbulence to stochastic models

The purpose of this section is to place this work in the continuity of a story which began around 1922 when Lewis Fry Richardson 10 described a turbulent flow as a cascade of vortices of all possible scales. Richardson summarized this description by a famous poem, see figure 2(a). This poem describes the energy transfer from the large scale (e.g., the size of a tea spoon when stirring tea) where energy is injected to the smallest scale where the energy is dissipated by viscosity in the form of heat: big whirls break down to smaller whirls, see figure 2 Why a physicist became so lyrical ? Certainly due to his difficulty in understanding and modeling turbulent flows: even today there is still no complete theoretical model of turbulent flows...

In fluid mechanics, the rate of turbulence of a flow is quantified by the Reynolds number Re. [START_REF] Frisch | Turbulence. The legacy of A. Kolmogorov[END_REF] This number Re is the ratio between the non-linear term and the viscosity term of the Navier-Stokes equations that govern the fluid dynamics. As Re increases, more and more instabilities appear and the flow becomes more and more disordered. For Re ≫ 1, the number of degrees of freedom is proportional to Re 9/4 (e.g., Re ∼ 10 10 ⇒ Re 9/4 ∼ 10 23 ). For such highly turbulent flows, the Navier-Stokes equations are difficult to solve analytically (because of serious mathematical difficulty † ) as well as numerically (because of too high numerical cost).

An alternative approach is to consider a turbulent flow as a disordered system to be modeled directly from experiment and data analysis. Such a model must of course obey known theoretical properties of turbulent flows as well. This approach is mainly based on the observation that although the detailed properties of a signal measured in a turbulent flow appear to be unpredictable, its statistical properties are reproducible. This observation motivated the search for a probabilistic description of turbulence as soon as in the 1930s. [START_REF] Taylor | Statistical theory of turbulence[END_REF]

Intermittency, multiplicative cascades and multifractals

In 1941, the famous probabilist A.N. Kolmogorov 12 proposed what is now referred to as the K41 theory. In summary, this theory consists of additional hypotheses to Navier-Stokes equations. We will focus on two of these hypotheses only. The first assumption is that all the possible symmetries of Navier-Stokes equations, usually broken by the mechanisms producing the turbulent flow, are restored in a statistical sense at small scales and away from boundaries in the limit of infinite Reylnolds number: in particular, homogeneity and isotropy are restored. The second assumption is that the flow is self-similar at small scales, i.e. it possesses a unique scaling exponent H for velocity differences δv(r, λℓ) law = λ H δv(r, ℓ). The velocity increments δv(r, λℓ) = v(r+ℓ)-v(r) must moreover obey the exact four-fifths law, derived from the exact Kármán-Howarth 13 (1938) relation, δv(r, ℓ) = -4/5εℓ where ε is the dissipation rate. Thus K41 theory predicts that H = 1/3. As a result, one may say that the K41 theory reminds to model the velocity field of a turbulent flow by a fractional Brownian motion with Hurst exponent H = 1/3. This model is described by a unique Hölder singularity exponent (see section 3.1) so that it is called monofractal. At this point it is interesting to note that Kolmogorov had just published in 1940 the article [START_REF] Kolmogorov | Wienersche spiralen und einige andere interessante kurven im Hilbertschen raum[END_REF] which originated the later work by Mandelbrot & van Ness 15 (1968) who defined the fractional Brownian motion (fBm) B H as presently used. Recall that the fBm B H (x) is the only self-similar Gaussian process with stationary increments. Its (generalized) Fourier spectrum is ∝ 1/k 1+2H ; for H = 1/3, one recovers the 1/k 5/3 spectrum usually observed in turbulence.

Meanwhile, in the 1940s, several works were devoted to the statistical description of fragmentation processes (e.g. Landau, [START_REF] Landau | On the energy loss of fast particles by ionization[END_REF] Ambarzumian 17 ) which share some similarity with the energy cascade from large to small scales: the phenomenology of Richardson's cascade evokes the fragmentation process of big whirls in smaller whirls. † The Clay Mathematics Institute of Cambridge, Massachusetts (CMI) has named seven Millenium Prize Problems: one of these problems is to make substantial progress toward a mathematical theory which will unlock the secrets hidden in the Navier-Stokes equations. Landau precisely made an objection to K41 theory by pointing out that it assumes the energy dissipation rate to be constant and uniform, although this assumption is not physical: the dissipation rate has to fluctuate in time and space. Landau made this objection in a footnote of the original edition of his book 18 published with Lifshitz in 1944. In brief, it means that the energy transfer is random and cannot be considered as uniform and constant, see figure 3(a). This is the origin of the so-called intermittency phenomenon.

In 1962, Kolmogorov [START_REF] Kolmogorov | A refinement hypothesis concerning the local structure of turbulence in a viscous fluid at high reynolds number[END_REF] (again !) and Obukhov [START_REF] Obukhov | Some specific features of atmospheric turbulence[END_REF] brought an enlighting answer to the problems raised by Landau. To this aim they introduced a log-normal dissipation rate ε ℓ at scale ℓ such that var [log ε ℓ ] = A + µ log(ℓ o /ℓ). This was the first model that took the intermittency phenomenon into account. It was also suggesting to build intermittency models based on the dissipation, see figure 3(a). This is the origin of the first multiplicative cascades introduced by the Russian school, e.g. Novikov, [START_REF] Novikov | The intermittency of turbulence and the spectrum of energy dissipation fluctuations[END_REF] Yaglom, [START_REF] Yaglom | Effect of fluctuations in energy dissipation rate on the form of turbulent characteristics in the inertial subrange[END_REF] . . . and later studied in finer details by Mandelbrot, 23 see figure 2(c). Mandelbrot remarked that the intermittency of multiplicative cascades could be described thanks to the notion of fractal dimension: in the limit of an infinite cascading process, the dissipation concentrates on fractal sets. Kahane & Peyrière 24 (1976) proved the mathematical theorems based on the theory of martingales supporting Mandelbrot's conjectures. These objects are described by a set of Hölder singularity exponents: they are multifractal.

From binomial to infinitely divisible cascades

The 1980s were rich in new results as far as multiresolution, multifractal, and intermittency notions are concerned. Wavelets 25 were "officially" born in 1984 , and pretty much at the same time multifractal analysis [START_REF] Hentschel | The infinite number of generalized dimensions of fractals and strange attractors[END_REF][START_REF] Frisch | Fully developed turbulence and intermittency[END_REF] was elaborated and the first accurate experimental observations of intermittency 28 were performed thanks to higher order statistics measurements. Therefore, the wavelet based multifractal analysis of the intermittency phenomenon in turbulent flows could be elaborated and developed, [START_REF] Jaffard | Multifractal formalism for functions, part 1 & 2[END_REF] in particular thanks to improved computer facilities for storage and analysis of huge amounts of data. [START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF][START_REF] Chainais | Intermittency and coherent structures in a turbulent flow: a wavelet analysis of joint pressure and velocity measurements[END_REF][START_REF] Arneodo | A wavelet-based for multifractal image analysis. i. methodology and test applications on isotropic and anisotropic random rough surfaces[END_REF] Section 3 recalls on wavelet based multifractal analysis.

In parallel with much data analysis about turbulent flows, 6 some people have also been searching for multifractal stochastic processes relevant for modeling. Multiplicative cascades are at the heart of this research. For instance, Schertzer & Lovejoy 33 (1987) introduced what they called "universal multifractals" build from multiplicative cascades using stable distributions to model rain and clouds. Arneodo et al. [START_REF] Arneodo | Random cascade on wavelet dyadic trees[END_REF] (1998) introduced multiplicative cascades on wavelet coefficients as an attempt to model signals such as the velocity signal in figure 3(b). We also want to mention the interesting purely mathematical work by Le Jan & Sznitman [START_REF] Jan | Stochastic cascades and 3-dimensional Navier-Stokes equations[END_REF] (1997) who proposed a solution in the weak sense to Navier-Sokes equations thanks to a construction using a multiplicative cascade.

Since 2000, various authors [START_REF] Schmitt | Stochastic equations generating continuous multiplicative cascades[END_REF][START_REF] Barral | Multiplicative products of cylindrical pulses[END_REF][START_REF] Muzy | Multifractal stationary random measures and multifractal random walks with log-infinitely divisible scaling laws[END_REF][START_REF] Bacry | Log-infinitely divisible multifractal processes[END_REF][START_REF] Chainais | Scale invariant infinitely divisible cascades[END_REF][START_REF] Chainais | On non scale invariant infinitely divisible cascades[END_REF][START_REF] Schmitt | A causal multifractal stochastic equation and its statistical properties[END_REF] have proposed a set of generalized versions of the Mandelbrot's discrete multiplicative cascades. These new constructions are called infinitely divisible cascades (IDC) because they rely on the use of infinitely divisible distributions. Section 4.5.1 recalls on the sub-family of compound Poisson cascades. From a mathematical point of view, these processes are defined as martingales that converge to the density of a multifractal (positive) measure. The main evolution compared to Mandelbrot binomial cascades is that such processes are truly stationary in time and share a true self-similarity with no special scale ratio (discrete cascades rely on some discrete underlying tree, e.g., the dyadix tree). Definitions have been extended to dimension D ≥ 1, 7, 43 see figure 6(b) for an example in dimension D = 2. It was recently shown that IDC were good candidates for the modeling of natural images. [START_REF] Chainais | Infinitely divisible cascades to model the statistics of natural images[END_REF] IDC therefore constitue a versatile framework for the modeling of measure-like physical quantities, e.g., turbulent dissipation, intensity of images... The fractional Brownian motion in multifractal time [START_REF] Mandelbrot | A multifractal walk down Wall Street[END_REF] then permits to model scalar quantities in one dimension with symmetric fluctuations such as one component of a turbulent velocity field as a function of one coordinate (e.g., v x (x), see figure 3(b)). The next step will be the definition of multifractal vectorial fields aiming at the modeling of a turbulent velocity field v(x, t).

The present work is the first attempt to model images of a turbulence like physical system, the quiet corona, by means of multifractal processes.

WAVELET MULTIFRACTAL ANALYSIS

Principle

Multifractal analysis aims at the characterization of the regularity of measures, functions or graphs of realizations of a stochastic process. [START_REF] Jaffard | Multifractal formalism for functions, part 1 & 2[END_REF] In the present context, the purpose is to quantify the presence of singularities in an image thanks to the multifractal spectrum D(h), where h is the so-called Hölder exponent that characterizes a singularity. In brief, f (x) is said to be locally Hölder regular with exponent h(x o ) at x o if h(x o ) is the highest exponent such that there exist a polynomial P (x) and a constant C with |f (x) -

P xo (x -x o )| ≤ C|x -x o | h for x in a neighborhood of x o . In nice cases, the polynomial P xo (x -x o ) is simply the Taylor polynomial of f (x) at x o . In the simplest case, P xo (x -x o ) = f (x o ) and the singularity is described by |f (x) -f (x o )| ∼ |x -x o | h .
The multifractal spectrum D(h) is defined as the fractal dimension of the set {x o |h(x o ) = h}. It can be defined for a measure or for a function.

Wavelets gather several interesting properties as a tool for multifractal analysis. [START_REF] Jaffard | Multifractal formalism for functions, part 1 & 2[END_REF][START_REF] Muzy | Multifractal formalism for fractal signals: The structure function approach versus the wavelet transform modulus-maxima method[END_REF][START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF] As a multiscale microscope, they precisely explore the local singular behavior of a function. If the function f has a singular behavior of exponent h in x o , then the wavelet coefficients T f (x o , a) located in x o at scale a are such that [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] |T f (x o , a)| ∼ a h as a → 0. This is indeed the case provided the wavelet has enough null moments to remain unresponsive to the local polynomial component P xo (xx o ) of f . First order increments f (x o + au) -f (x o ) would see the polynomial behavior only. In contrast, the regularity of a function may be precisely characterized thanks to wavelets even though a smooth polynomial behavior is superimposed to the singular part. Moreover, in its original formulation, the multifractal formalism applied to singular measures only. Wavelets have permitted the extension of the multifractal formalism to other objects like multifractal distributions.

In practice, the estimation of the local regularity of a multifractal function at every position x is numerically unstable. Therefore, an alternative approach is used which is based on the scaling behavior of global quantities called partition functions. In this approach, multifractal processes and multifractal images in particular are seen as a 'gas of singularities'. Multifractal analysis is then connected to some thermodynamical interpretation. In this view, the multifractal formalism is established when one can associate the multifractal spectrum D(h) to its Legendre transform ζ(q). Then q plays the role of the inverse of a temperature, D(h) is the entropy and ζ(q) is associated to the free energy. We will focus on this last function ζ(q) in particular. Indeed, the function ζ(q) is a set of multiscaling exponents of a partition function S(q, a) ∼ a ζ(q) usually built upon wavelet based averaged quantities (see next section). The exponents ζ(q) reflect the multiresolution statistics of a process. [START_REF] Castaing | Velocity probability density functions of high Reynolds number turbulence[END_REF][START_REF] Castaing | Fully developed turbulence : a unifying point of view[END_REF] For a stochastic multifractal process, they describe its higher order correlation structure. More precisely they contain information on the way the distributions of wavelet coefficients change from larger to smaller scales. This is the reason why we rather talk of 'multiscaling'. [START_REF] Chainais | Infinitely divisible cascades to model the statistics of natural images[END_REF][START_REF] Castaing | Fully developed turbulence : a unifying point of view[END_REF] The set of exponents ζ(q) will be used as a set of parameters of the models described in section 4.

Wavelet based estimates of the scaling exponents ζ(q)

Various methods based on wavelets are available to estimate the multiscaling exponents ζ(q). The Wavelet Transform Modulus Maxima (WTMM) method 32 is a priori very accurate for ideal multifractals but is quite delicate to carry out. Moreover, it needs large images to be accurate. This method consists of 4 steps:

1. Compute the 2D continuous wavelet transform of the image using derivatives of the Gaussian; 2. Extract chains of local wavelet modulus maxima; 3. Compute structure functions based on wavelet modulus maxima:

Z(q, a) = L∈L(a) sup (x,a ′ )∈L,a ′ ≤a M ψ [f ](x, a ′ ) q (1)
where ψ is the wavelet, L(a) denotes the set of maxima lines below scale a, L is a maxima line, and M ψ [f ] is a local maxima of the modulus of the wavelet transform of function f . 4. Estimate scaling exponents of these structure functions:

Z(q, a) ∼ a ζ(q) , a → 0. ( 2
)
When the power law scaling is obeyed over a sufficiently large range of scales, this method is very efficient. However, when applied to small images, it suffers from some sensitivity due to the step of the numerical determination of the maxima along maxima lines.

Recently another method was proposed that relies on similar ideas but applied to a discrete wavelet transform: Wavelet Leaders. [START_REF] Abry | Revisiting scaling, multifractal and multiplicative cascades with the wavelet leader lens[END_REF][START_REF] Lashermes | Wavelet leaders based multifractal analysis[END_REF] In this approach the supremum over maxima lines below scale a around x is replaced by the maximum over coefficients at octaves below j and within distance less than ±2 j of a given location k. Indexing for the 2D dyadic subsets λ ≡ λ (i) (j, k), i = 1, 2 or 3 correponding to the 3 possible orientations (horizontal, diagonal, vertical), j, k ∈ Z, the 2D wavelet ψ λ is essentially located around the subset λ. Let 3λ denote the subset λ and its adjacent dyadic subsets (at the same octave j). The wavelet leaders are defined as

L λ = sup λ ′ ⊂3λ |d λ ′ | (3)
where d λ ′ is the wavelet coefficient associated to the dyadic subset indexed by λ ′ . The modulus maxima are no longer taken over maxima lines as in the WTMM method but over sub-branches of the dyadic tree. The partition function becomes

S W L (q, j) = 1 N j Nj k=1 L q λ(j,k) ∼ 2 jζ(q) (4)
The more simple discrete wavelet transform (DWT) structure functions S DW T (q, j) are simply computed from the wavelet coefficients modulus, without taking any supremum:

S DW T (q, j) = 1 N j Nj k=1 |d(j, k)| q ∼ 2 jζ(q) (5)
where N j is the number of wavelet coefficients d(j, k) at octave j. This method is fast and efficient for positive values of q but becomes numerically unstable for q < 0 since the most probable value of the d(j, k) is zero.

Working with modulus maxima (WTMM or Wavelet Leaders) is theoretically the best approach that moreover gives the most stable numerical results. In particular, these methods permit to estimate scaling exponents ζ(q) for negative orders q. However, alike the WTMM method, the wavelet leaders method calls for a large scaling range [START_REF] Abry | Revisiting scaling, multifractal and multiplicative cascades with the wavelet leader lens[END_REF] which is not available in 256 × 256 images (only 8 octaves). Moreover, due to the use of a supremum over small scales, estimates remain sensitive to a possible cut-off of the scaling behavior at small scales. As a consequence, we have preferred to work with the more simple discrete wavelet based structure functions S DW T (q, j) so that the models proposed below will be fitted to experimental data on the basis of the ζ(q) for q ≥ 0 only.

STOCHASTIC MODELS FOR EUV QUIET SUN IMAGES

Toward multifractal models

We have selected a set of 54 images acquired at the wavelength 19.5 nm of EIT in 1997 dduring a period of minimal solar activity. We excluded images that clearly showed active regions, see figure 1(a), to keep only images featuring the quiet corona, see figure 1(b). Beyond classical analyses (histrograms, spectrum...) already studied in depth in the literature, these images exhibit a scale invariant property. [START_REF] Berghmans | Quiet Sun EUV transient brightenings and turbulence. A panoramic view by EIT on board SOHO[END_REF][START_REF] Delouille | Wavelet Spectrum Analysis Of Eit/Soho Images[END_REF] This can be illustrated for instance by the orientation averaged Fourier spectrum of Quiet Sun images, see figure 4(a). Moreover, the evolution of the distributions of wavelet coefficients through the scales is characteristic of a multifractal process, see figure 4(b). Multifractal analysis will permit to describe accurately the spatial correlation structure of these images. As a result, we identify two families of stochastic processes that reproduce a statistical behavior similar to Quiet sun images as well as a comparable visual aspect. Recall that we do not talk about direct and exact numerical simulations of the magneto-hydrodynamic flows in the Sun; we 'only' propose a statistical modeling of Quiet Sun images with appropriate stochastic models.

Experimental data

The images provided by EIT are of size 1024 × 1024 and represent the whole Sun, see figure 1. In order to limit the distortions due to the curvature of the Sun to about 15%, we have kept a 256 × 256 zone around the center of the apparent solar disk only. To avoid side effects, wavelet transforms have been computed on 512 × 512 images. To eliminate images that contained active regions, we have kept 512 × 512 images of which intensities I are such that P (I > U ) < 10 -4 . The threshold U was determined by computing the quantile U such that P (I > U ) = 10% for an image that contains an active region. This procedure yields a set of 54 images of size 512 × 512 without any strong active region in them and without any cosmic ray (usually seen as an isolated very brilliant point). The multifractal analysis is performed over the 512 × 512 images to avoid side effects in the wavelet transform but only coefficients associated to the central 256 × 256 square are taken into account, see figure 1(a).

Multifractal analysis of Quiet Sun images

Figure 5(a) shows the scaling behaviors of some structure functions estimated from the 54 Quiet Sun images. A linear behavior of log 2 S DW T (q, j) as a function of j = log 2 a is observed for 2 ≤ j ≤ 5. The exponents ζ(q) are estimated thanks to a linear regression is this scaling range.

Figure 5(b) shows the set of the resulting ζ(q) estimates for -1 ≤ q ≤ 5. Note that errorbars are computed as the empirical standard deviation of the set of 54 estimates. As remarked above, estimates for negative values of q are numerically unstable; this yields a large variability that results in wide errorbars. Estimates for q larger than 5 would not be reliable due to the set of data that is too small to guarantee sufficient statistics at higher obtained from the DWT structure functions (c) Exponents τ (q) deduced from ζ(q) by τ (q) = ζ(q)qζ(1).

orders. Errorbars are already quite important around q = 3, but it is also well known that our estimates might be biased when sample size is too small.

Figure 5(c) shows the exponents τ (q) = ζ(q) -qζ(1) of the measure underlying the images. A possible approach to model these images rely on the use of an integrated version of the density of some underlying multifractal measure described by a set of multifractal exponents τ (q) such that τ (0) = τ (1) = 0. Applying a 1/ k H filter modifies the multifractal exponents so that given τ (q) one gets a measure with multifractal exponents

ζ(q) = qH + τ (q). (6)
This is connected to the fact that the differentiation of a Hölder singularity of exponent h becomes a singularity of exponent h -1 while its integration yields a singularity of exponent h + 1. [START_REF] Mallat | Singularity detection and processing with wavelets[END_REF] It is then natural to consider the modeling of Quiet Sun images as a 1/ k H filtered version of a multifractal density. An interesting remark is that τ (1) = 0 by definition so that ζ(1) = H. This latter property will be useful to identify the parameters of the proposed models. The Fourier spectrum of a 2D multifractal density [START_REF] Chainais | Infinitely divisible cascades to model the statistics of natural images[END_REF] is ∼ 1/k 2+τ [START_REF] Schrijver | A new view of the solar outer atmosphere by the transition region and coronal explorer[END_REF] where τ (2) < 0 so that 2 + τ (2) < 2. The spectrum of the resulting model will be ∼ 1/k 2+τ (2)+2H which is consistent with the spectrum of Quiet Sun images that behaves as ∼ 1/k 2.9 .

Monofractal or multifractal ?

A basic modeling of Quiet Sun images would have taken into account the histogram of intensities and the Fourier spectrum. The exponential of a 2D fractional Brownian sheet (fBs) [START_REF] Bonami | Anisotropic analysis of some gaussian models[END_REF] B H appears as a good candidate to this purpose; H is called the Hurst exponent. The fBs can be approximately simulated in a very simple manner by applying a 1/ k H filter to a 2D Gaussian noise in the Fourier domain. The resulting process is Gaussian exhibits a 1/k 2+2H isotropic spectrum which is characteristic of scale invariance. Since the histogram of intensities in Quiet Sun images is concentrated on positive values and similar to a log-normal distribution, we may consider the exponential of the former 2D fBs. It has been usually observed when studying the statistics of natural images [START_REF] Chainais | Infinitely divisible cascades to model the statistics of natural images[END_REF][START_REF] Srivastava | On advances in statistical modeling of natural images[END_REF] that the power law behavior of the spectrum of a scale invariant process remains nearly unaffected by the exponentiation despite its non linear character. As a consequence, this model will exhibit both a 1/k 2+2H spectrum and a lognormal distribution of intensities. Estimates of the ζ(q) for the Quiet Sun prescribe H = 0.55 ± 0.06. One may be satisfied at first sight but this model is essentially monofractal. It has trivial multifractal exponents ζ(q) = qH, that is τ (q) = 0 in (6), in contrast with the observed ζ(q) which clearly exhibit a non linear behavior corresponding to a true multifractal behavior, see figure 5(b-c).

At this stage we still have to discriminate between two possibilities at least. As mentionned above, the first possibility is to use a 1/ k H filtered version of the density of a multifractal measure described by a set of multifractal exponents τ (q) which yields a process with multiscaling exponents ζ(q) = qH +τ (q). Several authors have introduced interesting definitions of such multifractal measures in 1 dimension. [START_REF] Schmitt | Stochastic equations generating continuous multiplicative cascades[END_REF][START_REF] Barral | Multiplicative products of cylindrical pulses[END_REF][START_REF] Muzy | Multifractal stationary random measures and multifractal random walks with log-infinitely divisible scaling laws[END_REF][START_REF] Bacry | Log-infinitely divisible multifractal processes[END_REF][START_REF] Chainais | Scale invariant infinitely divisible cascades[END_REF][START_REF] Chainais | On non scale invariant infinitely divisible cascades[END_REF][START_REF] Schmitt | A causal multifractal stochastic equation and its statistical properties[END_REF] These definitions have been recently extended [START_REF] Chainais | Multidimensional infinitely divisible cascades. application to the modelling of intermittency turbulence[END_REF][START_REF] Schmitt | On causal stochastic equations for log-stable multiplicative cascades[END_REF][START_REF] Chainais | Infinitely divisible cascades to model the statistics of natural images[END_REF] to dimension D ≥ 2 for image modeling purpose in particular. All these definitions are connected to the family of infinitely divisible cascades, a general framework for the multiplicative cascades originally introduced for the modeling of the dissipation field in turbulent fluid flows. [START_REF] Yaglom | Effect of fluctuations in energy dissipation rate on the form of turbulent characteristics in the inertial subrange[END_REF][START_REF] Mandelbrot | Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier[END_REF] Section 4.5 proposes a model according to this approach.

A second approach which is much less natural and intuitive but may be relevant as well cannot be completely discarded a priori. This is the special case where the multifractal behavior roots in the divergence of moments for orders q higher than some critical value. Linear fractional stable motions (in 1D) or sheets (in 2D) [START_REF] Samorodnitsky | Stable non-Gaussian random processes[END_REF] precisely exhibit such a special behavior. [START_REF] Chainais | Multifractal analysis and alpha-stable processes : a methodological contribution[END_REF] Thus the ζ(q) of these processes are non linear due to a phase transition at some critical order q = α: in this case ζ(q) is a piecewise linear function. Section 4.6 deals with this second approach to study its potential relevance.

A model based on compound Poisson cascades 4.5.1 Compound Poisson cascades

Compound Poisson cascades are a subset of the family of infinitely divisible cascades. They were originally introduced [START_REF] Barral | Multiplicative products of cylindrical pulses[END_REF] as a multifractal product of cylindrical pulses. This definition can be reformulated as a multiplicative cascade [START_REF] Chainais | Scale invariant infinitely divisible cascades[END_REF] and extended 7 to dimension D ≥ 2. The density Q ℓ (x) resulting from a compound Poisson cascade is defined by

Q ℓ (x) = i W f ( x-x i r i ) i IE i W f ( x-x i r i ) i (7)
where ℓ > 0 is some small limiting scale (a kind of resolution); the W i are independent identically distributed (i.i.d.) non negative random variables;

(x i , r i ) is a Poisson point process in R 2 × [ℓ, 1] with density dm(x, r) = (4/πr 3 )dxdr; f (x) = I 1 [-1/2,1/2] (x)
in the basic definition can be replaced by some compact supported non negative function ‡ The integration kernel f plays the role of some geometrical object in the image. It may also be useful to attenuate small scales discontinuities or to take into account some geometrical features of the images under study.

In the limit ℓ → 0, compound Poisson cascades are the density of a scale invariant multifractal measure characterized by a set of multiscaling exponents τ (q) = q(IEW i -1) + 1 -IEW q (at least within a certain range of q, see below). The design of some τ (q) function to the purpose of modeling reduces to the choice of the distribution of the multipliers W i .

An interesting property of the process Q ℓ (x) is that it can be interpreted as the intensity I(x) resulting from the scattering of a uniform light by a random superposition of transparent cylinders of sizes {r i } placed above positions {x i } and with i.i.d. random transparency W i . The centers x i of the cylinders are uniformly distributed on the plane; the radii r i are distributed by a scale invariant 1/r 3 law; the distribution of the transparencies W i is determined by the choice of the function τ (q) which is directly associated to their second generating function. The intensity of one pixel is therefore the product of the transparencies of the cylinders: this is a multiplicative cascade. This presentation points to the resemblance between CPC and other classical approaches in image modeling where elementary objects of random sizes are distributed in space following a Poisson point process. [START_REF] Srivastava | On advances in statistical modeling of natural images[END_REF] See 7, 44 for more details.

parameters H and τ (q)

The identification proceeds in 2 steps : first identifying H, and second evaluating τ (q) in equation (6). The parameter H describes the linear trends of ζ(q). The function τ (q) is a non linear concave function obeying τ (0) = τ (1) = 0 and controls the multifractal behavior of the final process. As a consequence, one has H = ζ(1) so that we will use for H the estimate of ζ(1) obtained from the multifractal analysis performed in section 4.3. Then τ (q) is estimated as τ (q) = ζ(q) -qζ (1). We get H = 0.55 ± 0.06.

We have studied the adequation of several compound Poisson cascade models to experimental data. The best fit seems to be obtained with a model such that W i = ((1 + T ) 1/T u) T where u is uniformly distributed in [0, 1]. Then one has: with T = 0.85 (see [START_REF] Chainais | Multidimensional infinitely divisible cascades. application to the modelling of intermittency turbulence[END_REF][START_REF] Chainais | Infinitely divisible cascades to model the statistics of natural images[END_REF] for a detailed presentation of available models). This theoretical expression is actually fine for q ≤ q * + in estimates based on the DWT structure functions S DW T (q, j). The existence of an upper bound q * + is due to the linearization effect. [START_REF] Lashermes | New insights into the estimation of scaling exponents[END_REF] In dimension 2, q * + is the value of q such that ζ(q) = 2 -qζ ′ (q). Therefore, one expects estimates of ζ(q) that behave as

τ CP C (q) = 1 - (1 + T ) q (1 + qT ) (for q ≤ q * + ) (8
ζ CP C (q) =    qH + 1 - (1 + T ) q (1 + qT ) for q ≤ q * + , 2 -qζ ′ (q * + ) for q ≥ q * + , (9)
where q * + ≃ 2.25 for T = 0.85. The estimates on figure 5(b) (black circles) are quite consistent with this theoretical prediction.

In practice the 1/ k H filtering is carried out in the Fourier space thanks to a fast Fourier transform and by using a 1/ k H frequential response that is truncated near the origin k = 0 since it is not defined at this point. We have chosen to impose a saturation at the value associated to the smallest available discretized frequency.

Once we have identified the parameters of this model we can numerically synthesize as many realizations as needed, at any desired resolution. To test the validity of this model, we compare the exponents ζ(q) and τ (q) estimated on 54 independent realizations of 512×512 model images to the exponents estimated in section 4.3, figure 5(b-c) on the 54 Quiet Sun images. Figure 5 shows that estimates on the images from the CPC model are nearly superimposed onto estimates from Quiet Sun images.

It is remarkable that the visual aspect of model images is quite similar to the visual aspect of Quiet Sun images, figure 6(b), except some geometrical structures which we did not aimed at featuring at this stage. We emphasize that only quantitative statistical parameters have been fitted to the results of the multifractal analysis and that no post-processing has been used. This approach is very systematic and theoretical, with no computer graphics artifacts. As a conclusion, this model seems to be quite relevant for EUV Quiet Sun images for λ = 19.5 nm. We recall moreover that such models receive some interesting physical interpretation, see previous section.

Remark: Figure 5 represents both the ζ(q) exponents and the τ (q) exponents (where the linear trend of the ζ(q) has been removed) to make comparisons more discriminating. Errorbars reflect the empirical standard deviation of the set of 54 estimates. They are not computed as an estimate of the variance of the estimator for one single image.

A model based on linear fractional stable sheets 4.6.1 Definition

In dimension D, the linear fractional stable sheet [START_REF] Samorodnitsky | Stable non-Gaussian random processes[END_REF][START_REF] Ayache | Local and asymptotic properties of linear fractional stable sheets[END_REF] (LFSS) of parameters α ∈ (0, 2] and H ∈ (H 1 , ..., H D ) ∈ (0, 1] D is the real valued symmetrical α-stable (SαS) process

L α = {L α (x), x ∈ R D } defined as L α (x) = R D D i=1 (x l -s l ) H l -1/α + -(-s l ) H l -1/α + dZ α (s), (10)
where {Z α (s), s ∈ R D } is the symmetrical α-stable Lévy sheet and x + = max(x, 0). The case α = 2 corresponds to the Gaussian case. Observe that, for every l = 1, ..., D, L α is a Linear Fractional Stable Motion (LFSM) in R of Hurst parameter H l along the direction of the lth axis. The choice H l = H, ∀l, ensures isotropy. The multifractal spectrum of LFSS is known to be very particular [START_REF] Chainais | Multifractal analysis and alpha-stable processes : a methodological contribution[END_REF][START_REF] Ayache | Local and asymptotic properties of linear fractional stable sheets[END_REF][START_REF] Jaffard | The multifractal nature of lévy processes[END_REF] since the multifractal properties of these processes are due to their high variability betrayed by the divergence of moments of order greater than α.

The scaling exponents ζ α,H (q) of an isotropic LFSS with parameter α and Hurst parameter H are given in D dimensions by: 58, 61, 62 which is a piecewise affine function. Lévy processes obey H = 1/α. Even though LFSS have no link with any underlying multiplicative cascades one can compute the corresponding τ (q) as for compound Poisson cascades as τ α (q) = ζ α,H (q) -qH which yields

ζ α,H (q) =      (H + D)q + D for q ≤ -1, Hq for -1 ≤ q ≤ α, (H -D/α) q + D for α ≤ q, (11)
τ α (q) =    Dq + D for q ≤ -1, 0 for -1 ≤ q ≤ α, -(D/α)q + D for α ≤ q. (12)
A remarkable feature of τ α (q) is that it is flat equal to zero for -1 ≤ q ≤ α which is never the case for any known multiplicative cascade. Let us remark that in dimension 2, the last linear function -(2/α)q + 2 is always below -q + 2 since 0 < α ≤ 2. Last, note that since we aim at modeling a positive quantity (light intensity) we used dissymetrical (β = 1) α-stable distributions in place of the usual symmetrical (β = 0) ones.

Identification of parameters H and α

Let us remark that the estimated exponents τ (q) are quite close to zero for -1 q 1.9. Moreover, they obey a quasi linear behavior for q ≤ -1 and q ≥ 1.9 both pointing to 2 at q = 0. This behavior is very similar to the multiscaling properties of the LFSS defined in previous section, see (12). Let us consider the adequation of some LFSS model to quiet Sun images.

Various wavelet based techniques 63 can be considered to estimate both H and α. Here the Hurst parameter H is estimated as H = ζ(1). Then α can be estimated in several ways. We have chosen to estimate α by best-fitting the function τ α (q) in ( 12) to the τ (q) estimated from Quiet Sun images. The estimated parameters that best fit the data are α = 1.9 and H = 0.55.

To validate this estimation, we compare the estimates from 54 independent realizations of LFSS similar to the Quiet Sun images. Simulations have been carried out by using an approximate Fourier method recently introduced by Scheffler & Biermé. [START_REF] Scheffler | Fourier series approximation of linear fractional stable motion[END_REF] In brief, the numerical methods consists of 2 steps: generating an α-stable Lévy image (i.i.d. random intensities) and Fourier filtering by 1/ k 2+H-2/α . Figure 5 shows the ζ α,H (q) and τ α (q) estimated from the LFSS model (triangles) compared to the ones estimated from Quiet Sun images (squares). The estimated ζ α,H (q) are within the errorbars of estimates on Quiet Sun images. This LFSS model might be relevant. although we get a less satisfying visual aspect compared to the previous compound Poisson cascade model (still no graphical post-processing has been made), see figure 6(c). The physical interpretation is much more delicate than that of compound Poisson cascades.

CONCLUSION

We have proposed a systematic approach to the statistical modeling of EUV Quiet Sun images from EIT onboard the SoHO mission (ESA/NASA). The purpose of this modeling is to be able to propose an extrapolation of what kind of images one should expect at resolutions much finer than what is actually available with present instruments. The main thrust of this work was to use the scale invariance property of quiet corona images. Therefore, we have first performed a wavelet based multifractal analysis of 54 Quiet Sun images at λ = 19.5 nm. Due to the small size of images (512x512), we have used Discrete Wavelet Transform partition functions rather than other more sophisticated tools like WTMM or Wavelet Leaders. As a result we get a set of multiscaling exponents ζ(q) and τ (q) that we further use as parameters to design stochastic processes similar to Quiet Sun images.

We have considered three families of models. The first model, the exponential of a fractional Brownian sheet, is a self-similar process built to simply reproduce the Fourier spectrum and the close to lognormal distribution of quiet Sun images. It would correspond to a monofractal behavior. The multifractal analysis unambiguously reveals the multifractal nature of Quiet Sun images. As a consequence we turned to multifractal models. The family of fractionnally integrated compound Poisson cascades belongs to the family of the more general infinitely divisible cascades. The parameters ζ CP C (q) can be nearly superimposed to those of Quiet Sun images. The visual aspect is moreover quite satisfactory, even though some geometrical features have not been taken into account. The third family of models is that of Linear Fractional Stable Sheets. Even though the parameters ζ α,H (q) are within errorbars of estimates from the quiet Sun images, they seem to fit the data in a worse manner and display a less satisfactory visual aspect.

If the α-stable was the best one, an important consequence would be that the images may have infinite variance. This should then be taken into account with care in usual image processing which usually assume a finite variance. However, even though we cannot definitely discriminate between the two multifractal models we have proposed here, it seems that the model based on compound Poisson cascades better fits our observations. At the moment, the accuracy of our numerical estimates is limited by the number of available data.

This study gives clear evidence that Quiet Sun images are highly intermittent and can be reasonably modeled by multifractal stochastic processes. This modeling will allow us to simulate images at any resolution by assuming the scale invariance property. Thus, we will be able to study the evolution of quantities such as the local signal to noise ratio and the local contrast as the resolution gets finer. This will help to prepare and calibrate future missions aiming at high and very high resolution imaging of the Sun.

Figure 1 .

 1 Figure 1. Examples of images from EIT at wavelength λ = 19.5 nm (a) with an active region, (b) without any active region, the Quiet Sun.

Figure 2 .

 2 Figure 2. (a) Poem by L.F. Richardson (1922); (b) corresponding phenomenological picture; (c) construction of a binomial multiplicative cascade.

Figure 3 .

 3 Figure 3. (a) Example of a turbulent energy dissipation rate ε(t) in a turbulent flow. (b) Example of a measurement of one component vx(t) of a turbulent velocity field.

Figure 4 .

 4 Figure 4. (a) Orientation averaged power spectrum for a given spatial frequency norm |k| from EIT 19.5 nm Quiet Sun images of 1997.(b) Normalized histograms of 2D discrete wavelet coefficients computed from 54 Quiet Sun images taken by EIT at 19.5 nm during the year 1997. (c) Histogram of intensities estimated from the 54 images.

Figure 5 .

 5 Figure 5. (a) Structure function (q = 0.33, 1, 2) estimated from DWT of Quiet Sun images in 1997. (b) Exponents ζ(q)obtained from the DWT structure functions (c) Exponents τ (q) deduced from ζ(q) by τ (q) = ζ(q)qζ(1).

  ) ‡ This may rejoin the random wavelet expansions evoked in Mumford & Gidas 59 (2001).

Figure 6 .

 6 Figure 6. Examples of 512×512 images of (a) Quiet Sun, (b) a fractionnally integrated compound Poisson cascade, (c) a linear fractional stable sheet (LFSS), (d) the exponential of a fractional Brownian sheet (H = 0.55).

ACKNOWLEDGMENTS

V. Delouille thanks the University Blaise Pascal of Clermont-Ferrand for the 1 month stay (may 2005) during which this work was initiated. This work has been supported by a France-Belgium Hubert Curien grant "Tournesol (Communauté francophone)".