A Century of Turbulent Cascades and the Emergence of Multifractal Operators

A century of cascades and three decades of multifractals have built up a truly interdisciplinary framework that has enabled a new approach and understanding of nonlinear phenomena, in particular, in geophysics. Nevertheless, there seems to be a profound gap between the potentials of multifractals and their actual use. For instance, it seems ironic that multifractals have been mostly restricted to scalar‐valued fields, whereas cascades were first invoked for the wind velocity. We argue that this requires to proceed to new developments of the multifractal formalism and to the emergence of multifractal operators. This paper therefore aims to first simplify the introduction to the most recent developments based on the analysis and generation of multifractal fields with the help of the group property of the responses of a nonlinear system to a scale change. The generators of the multifractal operators are introduced with the help of symmetries as simple and basic as orthogonal rotations and mirror symmetries. This leads in a rather straightforward manner to the large class of Gauss–Clifford and Lévy–Clifford generators that combine a number of seductive properties, including universal statistical and robust algebraic properties. At the same time, we obtain new results on the entanglement of spherical and hyperbolic geometries, as well as on the existence of finite statistics of these cascades.


Introduction
The concept of turbulent cascade is usually traced back to the celebrated Richardson's quatrain (Richardson, 1922): Big whirls have little whirls that feed on their velocity, and little whirls have smaller whirls and so on to viscosity. In fact, this quatrain was preceded by another paper (Richardson, 1920) that announced in a given way the concept of "quasi-equilibrium" introduced by Kolmogorov (1941), that is, how the kinetic energy is transferred from large scales to smaller scales down to the dissipation scale. Under various forms, the concept of cascades has been inspiring since that time for most theoretical developments in turbulence, see for instance Monin and Yaglom (1975), as well as for many other nonlinear phenomena, for example, climate (Lovejoy & Schertzer, 2013;Royer et al., 2008), astrophysics (Labini & Pietronero, 1996), and city dynamics (Dupuy, 2017;Murcio et al., 2015), just to cite a few examples.
It is remarkable that both Richardson and Kolmogorov focused on the (vector) velocity field u(x, t). However, the vector nature of the velocity field has been nevertheless forgotten due to the introduction of a scalar observable called the (second order) structure function defined as the ensemble average or mathematical expectation (denoted by E [ . ]) of the square of the velocity space increments: Earth and Space Science 10.1029/2019EA000608 S 2 (Δx, t) =E [(u(x + Δx, t) − u(x, t)) 2 ]. It has the interest to be Galilean invariant (e.g., independent of the mean velocity) and can be further reduced to a scalar function of a scalar variable (the separation distance |Δx|) with the help of the hypothesis of "local isotropy" (Kolmogorov, 1941). This (double) scalar reduction had been quite helpful to obtain the first quantitative law of the turbulent velocity field, the famous "Kolmogorov's 2/3 law" for the scaling of this structure function with respect to the separation distance |Δx|.
It is also significant that the breakthrough of multifractals (Benzi et al., 1984;Halsey et al., 1986;Parisi & Frisch, 1985;Schertzer & Lovejoy, 1984) largely resulted from revisiting explicit cascade models of the (scalar) energy flux density to smaller scales (x, t) (Mandelbrot, 1974;Novikov & Stewart, 1964;Yaglom, 1966), as well as the resulting, nontrivial scaling of high-order structure functions (Anselmet et al., 1984), where the square of the velocity increment is replaced by the qth power of its norm (S q (Δx, t) =E[|u(x + Δx, t) − u(x, t)| q ]), as somewhat anticipated by Kolmogorov (1962) and Oboukhov (1962). Again, the scalar reduction was helpful, in particular, in getting the physical understanding that intermittency corresponds to the fact that more and more turbulent regions of the fluid are not only smaller and smaller but have a smaller and smaller fractal dimension. The latter means that the activity of turbulence can be decomposed in a hierarchy of embedded fractals sets, hence the terminology of "multifractals," each set supporting a given level of turbulence activity. This loose definition of multifractal intermittency is generic in the sense that it applies to any kind of (scalar) activity (e.g., city growth) and could presumably apply to vector fields but in a not straightforward manner starting with the fact that their norms are not all equivalent.
As emphasized by Tchiguirinskaia (2017, 2018), there is a surprising gap between the roles of the respective dimensions of the domain and codomain of a multifractal field, respectively, the sets over which it is defined and into which it is valued (Bourbaki, 2004). Codomains have been rather limited to be one dimensional (1D), whereas the domain dimension could be arbitrarily large. This restriction has first prevented to address the key question of complex component interactions of vector fields, whereas this was somewhat done for dynamical systems mainly with 1-D domains but with larger codomain dimensions. This also explains the recurrent assumptions of mathematically convenient but unrealistic symmetries instead of studying the nontrivial symmetries corresponding to these interactions. In a more general manner, this restriction to 1-D codomain has considerably prevented the development and application of multifractals and kept barely open the Pandora box of multifractals despite their 35-year anniversary (Schertzer & Tchiguirinskaia, 2018). This paper builds upon these recent criticisms of scalar limitations and developments towards multifractal vectors, and introduces in fact multifractal operators.

Symmetries and Nonlinear Equations
The main argument in favor of the scaling of the velocity field u (x, t) is that the Navier-Stokes equations are invariant under dilation/contraction of the space T of any arbitrary scale ratio , respectively, smaller/greater than 1. Indeed, it suffices to rescale other quantities by appropriate powers of (operator T * in equation (1) to obtain the same equation only rescaled by a common prefactor: where (x, t) is the forcing, p(x, t)the pressure, (x, t) the fluid density, and its viscosity. The main differences between various approaches are that the viscosity is taken as 0 (Parisi & Frisch, 1985), rescaled to preserve a constant Reynolds number (Sedov, 1993) or rescaled like an eddy viscosity  as presently done in equation (1). Independent of these differences, the link between the scale symmetry of the generating equations and that of their solutions is not as direct as it may appear at first glance. Indeed, there is a divorce between isotropic rescaling of the space and that of the velocity vector: there is no reason that the velocity vector will remain collinear to its original direction. We therefore argue that it is indispensable to take into account other symmetries, such as rotational and mirror symmetries, whereas rescaling has been considered to apply only to rotational invariant observables such as structure functions. Decades of practice have thus kept scale symmetry separate from other symmetries. We develop the opposite and disruptive approach to consider the whole algebra of the generators of symmetries of a cascade. Schematic of a discrete (in scale) multiplicative cascade, where larger structures break into smaller ones in a self-similar manner, and their content (the energy flux density for turbulent cascades) is forwarded to the daughter structures via random multipliers. This 2-D representation can be understood as a 2-D cut of a 3-D process. Adapted from Schertzer and Lovejoy (1988).

Discrete Cascades
To move from a quatrain to concrete intermittency models was an important, historical step, although it was done with unnecessary constraints. In the caricatural but pedagogical case of "discrete (in scale) cascade models" (Mandelbrot, 1974;Yaglom, 1966), "eddies" are represented by adjacent and disjoint d-dimensional (hyper) cubes Δ i n of size L∕ n (n ∈ N, n = n 1 , i ∈ [[1, n ]] d ) (We use the semiclassical notation [[a, b]] for the interval of all integers between a and b included.) obtained by the hierarchical and iterative division of an initial cube Δ 1 0 of size L with a constant ratio of scales 1 (greater than 1 and usually equal to 2; see Figure 1). At the nth cascade step, there are therefore n subcubes Δ i n , and the energy flux density n is taken to be strictly homogeneous on each eddy Δ i n , that is, it is a step function.
where 1 Δ i n is the indicator function of the eddy Δ i n (= 1 if x ∈ Δ i n , and = 0 otherwise). The energy flux density n−1 at step n − 1 is multiplicatively distributed to subeddies of the nth step: with the help of the multiplicative increment n whose step function structure is as follows: generates in fact that of n (equation 2). Its components i n are usually assumed to be identically and independently distributed random scalars, as well as to be independent of the variables i n , that is, they are independent realizations of a given random variable ( d = denotes equality in distribution): 10.1029/2019EA000608 Despite their caricatural geometry, these models were helpful to get some key insights that we do not aim to fully review. Let us just mention that the convergence of the process relies on the "canonical conservation" of the energy flux, that is, a conservation on the average: n is said to be "conservative." It requires the necessary but not sufficient condition: Let us emphasize that despite the simplicity of the step functions involved at each cascade step, there is no pointwise convergence for n → ∞ but only a convergence to a measure that is furthermore singular with respect to the d-dimensional Lebesgue measure. This explains why we have to deal not only with pointwise functions but also with measures or generalized functions. Let us also point out that a major limitation of discrete cascades is that their tree structure imposes a particular metric: the 1 -adic ultrametric defined by the rank of the first distinct digit (in basis 1 ) between two eddy centers. This ultrametric generates a weird topology: two eddies of similar size can have only a far remote ancestor eddy, although being very close to each other for the usual (Euclidean) metric, but far away for the ultrametric. This shows that both metrics are incompatible. Consequently, translation invariance does not hold for discrete cascades, which is a serious drawback that is avoided by continuous (in scale) cascades (Schertzer & Lovejoy, 1987).

Cascade Generator, Singularities, and Continuous Cascades
The multiplicative nature of the discrete cascades invites to consider them as exponentials of additive processes: where Γ n is called the generator of the cascade that satisfies and in agreement with equation (3): where n has a structure similar to that of n , although obtained with the help of a logarithm that involves only indicator functions: where i n , similarly to the i n , are identically and independently distributed, that is, they are independent realizations of a given random variable : An important property of the cascade generator Γ n is that its first characteristic function (in the Laplace sense) Z Γ n (q) corresponds to the qth order statistical moment of the energy flux density n : where K Γ n (q) is the second characteristic functions (in the Laplace sense). These functions are also called, respectively, the moment and cumulant generating functions (still in the Laplace sense), whereas it is more usual to use the Fourier transform (that introduces an i prefactor in the right hand side of the first equality of equation 12). The mutual independence of the i n yields a logarithmic divergence of the generator Γ n with respect to the scale ratio n : and the scaling of the statistical moments of the energy flux density: 10.1029/2019EA000608 more precisely its multiscaling as soon as the scaling moment function K (q) is nonlinear. The i n are called singularities because they correspond to the rate of (algebraic) divergence of the energy flux density n with respect to the scale ratio n . There are stronger properties resulting from the facts that Z Γ n (q) is the Laplace transform of the probability distribution and that the power law is transformed into a power law by a Laplace transform or its inverse. These properties yield (see for instance Schertzer & Lovejoy, 2011) a scaling law for the probability distribution of the energy flux density n : where the function c( ) is a (statistical) codimension, as it could be understood with the help of a frequency interpretation of the probability: where the number of eddies having an energy flux density i n greater than n has a fractal dimension D( ), and d is the dimension of the embedding space.
With the help of these definitions, it is straightforward to establish that the duality of the statical moments and the probability distribution with respect to the Laplace transform yields a duality of moment orders and singularities with respect to the involutive Legendre function: This property was first underlined by Parisi and Frisch (1985) in the context of deterministic multifractals, a context that requires some ad hoc hypotheses that are not necessary for stochastic multifractals. When both c( ) and K(q) are differentiable, the previous equations can be rewritten like: which present a more explicit one-to-one correspondence between singularities and moment orders, as well as the fact that both curves c( ) and K(q) are envelopes of each other. A typical example of moment scaling function K(q) is that of conservative "universal multifractals" (Schertzer & Lovejoy, 1987;, which depends on the parameters C 1 ≥ 0 and ∈ [0, 2]: This expression will be derived in the more general context of Clifford operators (equations 60 and 94). The Legendre transform between K(q) and c( ) (equation 18) implies that: that is, C 1 is both the singularity and the codimension of the mean field. It therefore measures the mean intermittency of the field. On the contrary, being defined by the convexity of K(q) (of c( ) as well), and it therefore measures the variability of the intermittency with respect to the statistical order q or the singularity , hence its multifractality. In particular, the monofractal case of the so-called model corresponds to = 0, which is its lower bound (due to the convexity of K(q)), whereas = 2 is the upper bound corresponding to a Gaussian generator, as discussed in section 7.3.1.
Let us now emphasize that continuous values in scale cascades are based on the extension to continuous values of the scale ratio of two key properties of discrete cascades underlined above: • The energy flux density is the exponential of generator Γ .
• The generator Γ has a logarithmic divergence with the scale ratio .
Furthermore, these properties are sufficient to ensure the extension of most others, such as generalizations of the multiscaling of the statical moments (equation 14). Not surprisingly, multifractal operators will also be defined by both these properties or equivalently by generalisation of equation (14).
10.1029/2019EA000608 Figure 2. Commutative diagrams (see text) displaying (a, left) the pullback transform T * on the codomain of the field generated by the scale transform mapping the domain X; and (b, right) the push forward transform on the dual space of measures or generalized functions s generated by the pullback transform that maps the space of test functions. It is worthwhile to note that T * has the same role with respect to T * , as T has with respect to T * , although on different spaces.

Why to Leap from Scalars to Operators
This quick presentation of discrete cascades (section 3.1) points out that they are not a priori limited to (real) scalars, but they require a structure that has both addition and multiplication. Usual vector spaces have only multiplication by scalars and not by vectors. Therefore, there is no multiplicative cascade of ("usual") vectors. This may partly explain the lack of developments of multifractal vectors. A richer mathematical structure is required, that of algebra, that is by definition vector spaces possessing also a multiplication operation of vectors by vectors. The condition of convergence (equations 6 and 7) requires that such an algebra has a unity, that is, is a unitary algebra. Typical examples of such algebra are those of linear operators and of their matrix representations. These operators define cascades of vectors first by being vectors themselves. There is also the possibility to consider cuts by given subspaces of these operator cascades, as done in section 7 or to make them operate on a (usual) homogeneous vector field, as performed by Schertzer and Lovejoy (1995) for a discrete complex cascade applied to a strictly homogeneous 2-D flow.
However, we have also to note that discrete operator cascades face not only similar limitations to those of discrete scalar cascades but also some new limitations. Indeed, products of one-parameter operators in general do not provide a one-parameter group. Furthermore, they may not preserve all their degrees of freedom. This is the case for the representation of 3-D rotations with the help of Euler angles, where the so-called gimbal lock phenomenon can reduce the degrees of freedom from 3 to 2.
To overcome these limitations, section 3.2 points out that we have to consider exponentials of operators rather than their multiplications, in analogy with how scalar-valued continuous cascades were obtained. Before doing this, it is worth to pause and to quickly evaluate what can bring this leap from scalars to operators. A particular novelty, similar to that at the emergence of quantum physics, is that operators in general do not commute, contrary to scalars. There is therefore no surprise that noncommutative algebra will occupy a central place in the sections that follow, particularly section 6, but we have first to define the response of a system to a change of scale.

Responses to a Scale Change
In this section, we consider the responses of a field and a (mathematical) measure or generalized function to a change of scale T . Both responses are needed due to the expected singular limit of cascades (section 3.1).

Response of Fields to a Scale Change: The Pullback Transform
Let be a field that maps the domain X (e.g., space-time) into the codomainX (e.g., velocity components), both being either vector spaces or manifolds: and a scale change operator T mapping the domain X into itself. In the simplest case, it corresponds to a contraction/dilation for any given positive scale ratio ( > 1 for a contraction and 0 < < 1 for a dilation): to the push forward transform T * , , by duality. It is worthwhile to note that the pullback transform T * has, again, a similar role on C(X,X) as the scale transform T on X.
We will discuss in section 5 strongly anisotropic extensions preserving the multiplicative group property with respect to the scale ratio : Independent of its precise definition, the map T of X into itself defines a map T * ofX into itself in agreement with: This equation, as illustrated by Figure 2a, means there are two paths to obtain the same result, a property called "a commutative diagram." The response T * to T is called the "pullback transform" (or "composition operator"; Shapiro, 1993): it pulls back the map from the new coordinates = T x to the old coordinates x, with respect to the "change of coordinates" T , hence a generalization of the contravariant change of coordinates. The star superscript is indeed in agreement with the usual notation of contravariant tensors.

Examples of Pullback Transforms
The simplest case corresponds to the scale invariance of a (geometric) fractal set A of indicator function 1 A (Id denotes the identity transform): One may note that if A has a finite outer scale (e.g., A is compact), one must consider the intersections with a small enough sphere S to eliminate finite side effects. The classical definition of a fractal process B(t) with stationary increments (Lamperti, 1962) is rather similar to a structure function but it makes intervene an equality in distribution therefore for all moments: which corresponds to: The scaling exponent H refers to the hydrologist Hurst (Hurst, 1951;Klemes, 1974;Koutsoyiannis, 2002), who uncovered this exponent from empirical analyses of long Nile flood series. The value H = 1∕2 corresponds to the Brownian motion and normal diffusion, whereas other values correspond to fractional Brownian motions (Kahane, 1997;Mandelbrot & Van Ness, 1968) and other anomalous diffusion processes (Lévy, 1965;Painter, 1996;Yanovsky et al., 2000). The (relative) generality of these additive processes unfortunately gave credence to the uniqueness of the scaling exponent.

Response of Measures or Generalized Functions to a Scale Change: The Push Forward Transform
By duality, the pullback transform generates the push forward transform T * , (see Figure 2b for illustration), which maps into itself the set of (mathematical) measures or generalized functions , that is, the dual space C(X,X) ′ of the linear transforms of a given space C(X,X) of continuous test functions mapping X intoX: which again corresponds to a commutative diagram (Figure 2b). It is worthwhile to note that the pullback transform T * acts on C(X,X) as the scale transform T does on X: Due to the strong similarity between T * and T * , , we will use the common notationT to represent both of them when their differences are not relevant. The duality implies that the pullback transform T * indeed pushes forward the measure from the old coordinate x to the new coordinate T x and thus generalizes the covariant coordinate transform of (finite dimensional) dual vector spaces. A straightforward example is that of a fractal measure A , that is, a generalization of the Dirac measure such that it has unit mass homogeneously distributed on a fractal set A of dimension D A : It is worthwhile to note that the renormalization group approach (Wilson, 1971) and its variants define transformations similar toT by a decimation process (i.e., reducing the degrees of freedom of a system) rather than by a (geometric) change of scale.

Generalized Scale Transform and Metrics
We recall here the possibility and necessity to generalize the concept of scale transform T well beyond isotropic dilations/contractions. This will be done by (i) preserving the multiplicative group property of T (equation 24); and (ii) defining an associated, generalized notion of scale. Before discussing the latter, let us note that the group property of T propagates in a straightforward manner to the pullback transform T * , and then to the push forward transform (see Figure 3 for an illustration).
Although the fractal geometry (Mandelbrot, 1977(Mandelbrot, , 1983 claims to be fundamentally non-Euclidean, it relies on the Euclidean metric, which is isotropic. However, there is an abundant evidence of anisotropic fields and patterns, notably in geophysics, that exhibit an extreme variability over a wide range of scales. It is therefore indispensable to avoid the usual hypothesis of rotational symmetry prior to a scale symmetry (Schertzer & Lovejoy, 1985b, 2011. Anisotropic cascades and metrics therefore have become indispensable. For instance, considering only (quasi-isotropic) 2D and 3D regimes for atmospheric turbulence has been dramatically misleading for decades because none of these regimes seem to be relevant, contrary to a unique regime of Figure 5. (a, left) I is the orthogonal rotation (with respect to the origin of axes), and (b, right) J and K are the axial/mirror symmetries with respect to the first bissectrix and abscissa axis, respectively. It is easy to check that where 1 denotes the identity, as well as J = IK = −KI and IJK = −1 and therefore equation (40). (Schertzer & Lovejoy, 1983, 1985a, 1985bSchertzer et al., 2012). The corresponding anisotropic metric ||.|| is illustrated by Figure 4 that displays balls B defined by it: whereas the corresponding balls defined by the Euclidean metric | · | will be self-similar spheres. Despite this strong difference, both metrics have the common properties of nondegeneracy and nonnegativeness: and also homogeneity, although no longer limited to isotropic changes of scale: However, the triangular inequality (or subadditivity property) is no longer required for the generalized scale || · || but only the weaker property of ball embedding (Schertzer et al., 1999): The simplest case of generalized scale invariance corresponds to that of linear GSI, where the scale transform is a multiplicative group generated by a given matrix G (= Id, for isotropic scaling): Diagonalizable generators G with a positive spectrum Spec(G) = { i } and corresponding eigenvectors e i yield straightforward examples of generalized metrics for any p > 0, generalizing the classical L p norms, which correspond to G = Id (Spec(G) = {1}) and p ≥ 1 (For p ∈ [0, 1), the scale defined by equation 38 does not fulfil the triangular inequality even for G = Id): which obviously satisfies the homogeneity property (equation 35).
Equation (38) illustrates that, like for classical norms, generalized scales can correspond to the same scale transform, without being necessarily equivalent. The notion of equivalence between two generalized scales (|| · || 1 and || · || 2 ) remains the same as for classical norms:

Rotations, Mirror Symmetries, and Pseudo-Quaternions
Until now, we have focused on the scale symmetry. However, the previous section shows that the generator G of the scale transform T is not necessarily reduced to the identity with important consequences as already shown by the explicit example given above of generalized metrics (equation 38), although limited to a diagonal generator. Nontrivial generators G can be understood as resulting from a combination of symmetry generators. Indeed, such generators can be added and composed, generating therefore algebras, where the multiplication corresponds to the composition. Three simple, linear plane symmetries (I, J, and K) and their iterated applications (I 2 , J 2 , and K 2 ) are displayed in Figure 5: I is the orthogonal rotation (with respect to the axis origin), and J and K are the axial/mirror symmetries with respect to the first bissectrix and abscissa axis, respectively. These properties are summarized by: where 1 denotes the identity application, I is therefore a square root of the negative identity (−1), whereas J and K are square roots of the positive identity (+1). The action of I corresponds to a multiplication by the imaginary number i in the complex plane, and K corresponds to the complex conjugation. Both preserve angles, but only the former preserves their orientation, a property required by the strict definition of conformal transforms, whereas a mirror symmetry does inverse the angle orientation. The former transform is holomorphic, contrary to the latter.
Another profound difference between I and K is that I obviously results from the iteration of a small (possibly infinitesimal) rotation, that is, the tip of the vector follows from the position X 0 to X 1 a quart of a circle, which is a geodesic defined by I. A similar property is not so obvious for K: indeed the corresponding geodesic is a hyperbola and the corresponding rotation is no longer a classical, spherical rotation but a less familiar, hyperbolic rotation, as discussed in section 6.4. These properties also hold for J because, as it can be seen in Figure 5, the transforms I, J, and K are not fully independent: the composition of two of them yields the third one up to a ± sign. In particular, J and K correspond to each other via an orthogonal rotation (±I). Furthermore, due to this sign dependency on the order of the composition of these operators, they all anticommute: where {., .} and [., .], respectively, denote the anticommutator and the commutator: The commutator is a particular case of the celebrated Lie bracket of a Lie algebra. The previously highlighted properties of (I, J, K) show that: • Two among these three operators are sufficient to generate (by composition) the basis (1, I, J, K) of a larger vector space (and algebra). For reasons discussed below, we have often called this algebra quasi-quaternions (like Yaglom, 1968 andKocik, 2007) or pseudo-quaternions, although this term was also used by Okubo (1978) for a more involved nonassociative algebra. They have also been named coquaternions (Cockle, 1849) and more recently split quaternions (Rosenfeld, 1988). • Further compositions of these operators do not yield any new operator, that is, it is a closed algebra. • Due to its bilinearity and symmetry, the anticommutator defines a scalar product for which (I, J, K) are orthogonal due to equation (41) (with the notation (I, J, K) = (e 1 , e 2 , e 3 )): The equation that follows was slightly modified to be more accurate: • Which in turn defines the nonpositive definite quadratic form Q(v) =< v, v > for all vectors v's generated by (I, J, K). • A matrix representation of the symmetries (1, I, J, K) is: but this is by no means the unique representation of pseudo-quaternions, including linear ones.

Quaternions
Equation (40) that summarizes the pseudo-quaternion properties is very similar to the celebrated quaternion equation (Hamilton, 1844): This explains the aforementioned terminology. However, there is an important difference: I 2 , J 2 , and K 2 are all three square roots of the negative unity (−1), none of them of the positive unity (+1). Quaternions generate only spherical rotations. Three independent square roots of the negative unity cannot be obtained neither in the 2-D (real) plane, nor in the 3-D (real) space. In fact, it requires a 4-D (real) hyperspace so that one can combine two 2-D rotations or two 2-D mirror symmetries to obtain new squares of minus unity. This can be easily seen on the classical matrix representation of the quaternions expressed here as 2 × 2 block matrices of the pseudo-quaternions (equation 45):

Clifford Algebra
We introduced the algebra of pseudo-quaternions H ′ and of quaternions H (classical notations [The H symbol for the quaternion algebra should not be confused with the Hurst exponent H discussed in section 4.2.]) because they are both pedagogic examples of a more general type of algebra that have strong structural properties that we would like to use to generate symmetry groups for multifractal vector fields, in short to generate multifractal operators. H ′ and H are indeed both Clifford algebra (Baylis, 2004;Hagen & Scheuermann, 2001;Trautman & Warszawski, 2006), but their quadratic form, a basic feature of Clifford algebra, has a different signature as can be inferred from equations (40) and (46) and as discussed below. Indeed, a real Clifford algebra Cl p,q (R) is generated by n operators e i that are orthogonal (with respect to the anticommutator, equation 44), whose p square to the positive identity (+1) and q to its negative counterpart (−1). (The symbol q is classical for Clifford algebra but will be used only in this section in this sense to avoid any confusion with the [also] classical use of the same symbol for the statistical moment order in other sections.) They define therefore a quadratic form of signature (p, q), with p + q = n, for any operator generated by the operators e i such that: For instance, the algebra of the pseudo-quaternions H ′ is isomorphic to the linear algebra l(2, R) of 2 × 2 real matrices spanned by 1, I, J, and K (equation 45) and can be generated by V = {I, J} or {K, I}. Therefore, it can be denoted by Cl 2,0 (R), as well as by Cl 1,1 (R) because it can also be generated by {J, K}. On the contrary, the algebra of the quaternions H, spanned by 1 2 , I 2 , J 2 , and K 2 (equation 47), univocally corresponds to Cl 0,2 (R), although it can be generated by V = {I 2 , J 2 }, {J 2 , K 2 } or {K 2 , I 2 }, but all these spaces have the same signature (0, 2). One may note the simpler examples: Cl 0,0 (R) is isomorphic to R (V = ∅, no vector, and only scalars) and Cl 0,1 (R) to C (a unique generating vector I, which squares to −1, V = {I}), Cl 1,0 (R) seem to be nonclassical (generated by V={J} or {K}). Based on its generation, it is straightforward to demonstrate that the dimension of the Cl p,q (R) algebra is : It will be useful in what follows to consider the extension Q v of the initial quadratic form Q over the vector part of the whole algebra Cl p,q (R): and its signature (p v , q v ) (n v = 2 n − 1 due to equation 49), which is univocal contrary to the signature (p, q) of the generating algebra, as discussed above. For instance, (p v , q v ) = (2, 1) for H ′ (with (p, q) = (2, 0) or (1, 1)), and (p v , q v ) = (0, 3) for H (with (p, q) = (0, 2)). It is sometimes useful to also consider the extension of this quadratic form on the "scalar part" generated by the unity operator, that is, Q v (v 0 ) = v 2 0 . For simplicity's sake (as well as to avoid confusion on the meaning of the symbol q), a Clifford algebra will be merely denoted by the symbol .

Exponentials, Spherical, and Hyperbolic Rotations
We now discuss in more detail the questions of (infinitesimal) rotations and geodesics. It is convenient to introduce the terminology of unitary operator for square roots of either positive or negative identity, that is, u 2 = ±1, to investigate both the spherical and hyperbolic cases. We are primarily interested to see how is the transformed geodesic defined by u in the algebra (i.e., a straight line) by exponentiation into the corresponding group. Although, this transform is always a geodesic, it will strongly depends on the sign of u 2 . Indeed, one easily obtains, respectively, both are generalizations of the Euler-Moivre identity in the complex plane (with u = ), and the angle remains a curvilinear coordinate along a geodesic but in a hyperbolic geometry framework (Milnor, 1982) for u 2 = 1, instead of the (usual) spherical geometry for u 2 = −1. Like in relativity, both geometries are separated by the light cone defined by the nilpotent operators (also called isotropic) that square to 0 (u 2 = 0). The exponential series expansion is then limited to its two first terms:

Overview
Once a generator algebra  is chosen (in fact, a Clifford algebra), we can transform it into a stochastic algebra  p by defining generators from a given probability space onto . We do it by first considering the Gaussian case (section 7.2), that is, the Gauss-Clifford generators that have not only their own interest but they greatly facilitate the introduction of concepts necessary for the Lévy-Clifford generators (section 7.3, in particular Figure 7). The latter are very broad generalization of the former but require to overcome some technical difficulties that are already present for the (scalar) Lévy variables (section 7.3.1). Mainly, theoretical statistics are not always finite, and the corresponding empirical estimates are then sample dependent. We physically relate these statistical divergences to the underlying distribution of jumps that become more complex with 10.1029/2019EA000608 increasing codomain dimensions (section 7.3.2), as well as with the nonpositive definiteness of the quadratic form Q v (equation 50) of  and its associated scalar product < . > (see also equation 44): We precisely define admissible codomains of jumps  ↑ and corresponding domains of finite statistic orders  > , both having a cone structure. This enables to easily obtain the scaling moment function of Lévy-Clifford operators (section 7.3.3) and to conclude with examples, including with the help of a nonstandard generalization of Lévy vectors (section 7.3.4).

Gauss-Clifford Generators and Directional Statistics
A basic property of Gaussian vectors is that not only all their components are Gaussian variables, but the converse is also true: any random vector having all components to be Gaussian is a Gaussian vector. We will see a broad generalization of this property in section 7.3.2. This underlines the importance of the scalar product of the Clifford algebra, which defines the component Γ u of any (stochastic) Γ Clifford generator along any deterministic u operator: in particular, along any unitary vector u (u 2 = ±1) such as the vectors e i of the canonical basis  (equation 50) with Whereas statistics of the responseT = exp(Γ ) seem at first glance to be very complex, the components Γ ,u 's enable us to define scalar-like statistics ofT along any direction u that belongs to the same Clifford algebra. This enables us to proceed to a kind of polar decomposition. Indeed, we can define directional statistical moments of the operatorT for any order q ∈ R, as being the statistical moments of the component Γ ,qu = qΓ ,u : where the arrow ⃗ u denotes along the direction u. Similarly to the scalar case (equation 12), the last equality of equation (57) implies that the directional statistical moments of the responseT (the exponential of the generator Γ ) correspond to the first characteristic function (in the Laplace sense) Z Γ ,u (q) of the component of generator Γ ,u along this direction u due to: It is important to note that this shows, like for the scalar case (equations 13 and 14), that multiscaling of the directional moments E[T q⃗ u ] is obtained as soon as the second characteristic function K Γ ,u (q) has a logarithmic divergence with the scale ratio . This confirms that it could be taken as a definition of multifractal operators, as suggested in section 3.2.
Although much less used than the Fourier transform, the Laplace transform of a Gaussian variable X u of mean X u and variance 2 X u is straightforward to derive and yields the corresponding second characteristic function: This enables to define a Gaussian white noise "density" by its components u 's of variance 2C 1 (u) ≥ 0 in the direction u: The scale prefactor C 1 (u) of the Gaussian white noise density generalizes the codimension C 1 of the main field (equation 20) but with a (nonsurprising) vector direction u dependency. The integration of this white noise density over pixels of size , dimension d, and centers x i yields a white noise u : This "subgenerator" is then fractionally integrated to obtain a generator with a resolution logarithm divergence (equation 13): where the domain D( ) is limited by the external and internal scales, that is, denoting the ball of scale L, dimension d, and centered on x. One may note that B L (x) is in general not a sphere but a hyperboloid of revolution (like in Figure 6) bounded by its spherical basis of size L. G(x) is a nonnegative, scaling Green's function with parameters e and to be determined to obtain a scaling/group property This fractional integration can be discretized in the following manner: Due to the fact that variances of independent variables just add, Considering the limit → 0, we obtain (In fact, we are mimicking the classical Ito calculus.):

10.1029/2019EA000608
where B 1 is the d − 1 dimensional surface of the unit ball. Selecting the following prefactor e and exponent − of the Green's function (equation 63): and subtracting the centering term C 1 (u) log( ) yields: which shows thatT is conservative, that is, its average is the identity: which generalizes equation (6). As expected, the scaling moment function K Γ ,u (q) differs from that of a conservative, scalar multifractal with a Gaussian generator (equation 19) only by the direction dependency of the codimension C 1 (u) of the mean intermittency (equation 20).

Lévy -Clifford Generators 7.3.1. Lévy Stable Variables
We have basically used in the previous section the stability property of Gaussian variables that can be expressed in two slightly different manners that are easily shown to be equivalent: Indeed, if X is a Gaussian variable with mean and variance 2 (X ∼  ( , 2 )), they correspond, respectively, to: This stability is in fact equivalent to the central limit theorem that is, Gaussian variables are attractive for renormalized sums of independent realizations Y (i) (i = 1, n) of a random variable Y having a finite variance 2 (therefore a finite mean ). But, this theorem has been broadly generalized to Pareto variables Y (i) having infinite variance like X; more precisely, their common property is that they have a power law tail exponent ∈ (0, 2): which is the critical exponent of the statistical moment divergence.
Indeed, Gaussian variables are only a very particular case of the large class of the random variables called "Lévy stable variables" (Duan, 2015;Gnedenko, 1954;Feller, 1971;Kahane, 1985;Lévy, 1925Lévy, , 1965Zolotarev, 1986), which are defined by the equivalent properties of stability (equations 70 and 71) and attractivity (equation 73). Their main parameter is the "Lévy stability index" ∈ (0, 2] (2 for the Gaussian case), which is the exponent of the stability equation (equation (70)) and the inverse of the generator of the multiplicative rescaling group c(n) of equations (71) and (73): The exceptionality of the Gaussian case is that it has a fall-off faster than an exponential, and all its statistical moments therefore converge. Both behaviors correspond to = ∞ for equation (74), whereas = 2 is required for the renormalization (equations (71-73)). On the contrary, the divergence of moments of the Lévy stable variables has drastic consequences on their characteristic functions because of the trivial inequalities for any nonnegative random variable X: which shows that the characteristic function Z X (q) is infinite for any statistical order q > 0 for any random variable X with a finite power law tail exponent for its nonnegative values (equation 74 for max(X, 0) instead of |X|). This is due to the presence of jumps Y (i) whose intensity is Pareto distributed (i.e., has a power law distribution with an exponent ∈ (0, 2)). Only asymmetrical Lévy stable laws can therefore be used (Brax & Pechanski, 1991;Fan, 1989;Schertzer & Lovejoy, 1987, 1997, that is, those with only a power law tail for negative values that do not introduce any positive order moment divergence. The role of the jumps Y (i) are better understood with the help of a random version of the generalized central limit theorem (equation 73) that corresponds to a Poisson sum of the jumps Y (i) instead of a deterministic sum (and considering the limit of an infinite sample): where N is the number of jumps with a Poisson distribution of parameter (its mean, as well as its variance and all other cumulants), that is, with the probability and characteristic functions as follows: It yields in a rather straightforward manner: with the choice d(N) = (to be justified below). Having only negative jumps, Y (i) (0 < < 2 and ′ = − ): The last equality is obtained with the help of two successive integrations by parts and the choice d(N) = that corresponds to a regulation of e −q by the two first terms of its expansion to prevent the divergence of the integral at = 0 for < 2. It also establishes the relation between the average density of jumps and the corresponding codimension C 1 of the mean intermittency (Γ E denotes the Euler gamma function): both being nonnegative for ∈ (0, 2) (see Figure 8). Γ E (− ) has simple poles for nonnegative integer values of and yields: The behavior of ( ) for nearby 2 just confirms that the density of jumps smoothly vanishes for → 2. One may note that before Schertzer and Lovejoy (1987), the Laplace characteristic function has been scarcely used for Lévy stable variables and almost exclusively for the case < 1, where the probability is one sided, that is, corresponds to stable Lévy stable variables having an upper or lower bound (Feller, 1971).

Lévy Stable Vectors
Being linear, the properties of stability (equations 70 and 71) and attractivity (equations 73 and 77) can be immediately extended to vectors, the centering terms d and d(n) becoming vectors, whereas the renormalizing coefficients a, b, c, and c(n) are classically kept scalar (Lévy, 1937;Nikias & Shao, 1995;Paulauskas, 1976;Samorodnitsky & Taqqu, 1994). A more general and more involved case will be discussed in section 7.3.4. Still, due to this linearity, these properties hold for all components X u . Somewhat surprisingly, the converse is also true for "strictly" stable vectors, that is, having no centering terms d in equation (70) and d(n) in equations (71) and (73). This is due to the (strict) linearity of stability (equation 70 with d = 0), which applies both to the vectors X (1) , X (2) , and X and their components X (1) u , X (2) u , and X u , as well as the equality between their respective characteristic functions (either Laplace or Fourier, the latter adds only an i prefactor to the scalar product): Furthermore, the equality of characteristic functions is equivalent to equality in distribution.
Because the fundamental direction vector is qu in equation (83), the characteristic functions Z X u (q) and K X u (q) are invariant under the change of sign of both q and u: The required asymmetry of the Lévy stable laws to have finite statistics for positive orders (equation 76) together with this symmetry introduces a difficulty that does not exist for Gaussian laws, which are on the contrary fully symmetric. For Lévy stable variables, this asymmetry requires that the codomain of jumps Y 's is restricted to be R − . We now determine in any arbitrary dimension this codomain that we denote  ↑ and its opposite  ↓ These notations correspond to u ∈  ↑ points towards jumps, whereas u ∈  ↓ points in the opposite direction. Contrary to the 1-D case, the domain of positive orders yielding finite statistics  > , to be defined below, is not identical to  ↓ (=R + for the 1-D  = R). Obviously,  ↑ and  ↓ are asymmetric, linear cones (Recall that a real linear cone  is a set invariant by multiplication by any nonnegative real number, that is,  = R + ): where  ↑ 1 and  ↓ 1 (= − 1 ↑) are "unitary" sections of these cones (what is really required is that these sections are bounded). For instance, for  = R:  ↑ = R − and  ↓ = R + ,  ↑ 1 = {−1},  ↓ 1 = {+1}, but other real numbers of the same sign would be as convenient. Obviously, these properties are rather unchanged if a 1-D  ↑ is immersed into a larger space than R.
The projection Σ(N) u on a direction u of the Poisson sum Σ(N) of jumps Y 's integrates all the projections of radial (1-D) sums of these jumps (equation 80) along directions v's, which are distributed according to a measure that is discrete or continuous. Σ(N) and Σ(N) u have therefore the second characteristic functions as follows: with a straightforward generalization of equation (81) to the relation of the jump density (dv) and the codimension C 1 of the mean intermittency: 10.1029/2019EA000608 However, K Σ(N) u (q) is finite for q > 0 only if < u, v > is nonnegative in equation (87), that is, only if the direction u belongs to:  > being a cone, because both A and  ↑ are, we can often replace it by one of its sections  > 1 . For instance, in the 1-D example of  ↓ = R + {e} (i.e.,  ↓ 1 = {e}),  > is as follows: where (e) is the "half space" of  containing the vector Q v (e)e and limited by the "hyperplane" e ⟂ orthogonal to e: and where a section  > 1 of the cone  > is One may note that any isotropic/idempotent vectors e belong to their orthogonal set e ⟂ and therefore require a refined analysis.

Lévy-Clifford Generators and Operators
The characteristic function of the projection on a direction u of the Poisson sum of all jumps Y (equation 87) shows that a Lévy stable noise density generalizes the Gaussian white noise "density" (equation 60) is as follows: The corresponding subgenerator is readily obtained from equation (94): which is fractionally integrated, like in the Gaussian case (equation 62), with the help of a scaling Green's function(equation (63)). Systematically replacing the exponent 2 by in equations (66) and (67) yields the generalization of the directional moments of the conservative Gaussian-Clifford operatorT (equation 68) into that of Lévy-Clifford operators: after having introduced a normalizing factor exp(− C 1 (u) −1 log( ))) to impose the canonical conservation (equation 69).

Examples of Lévy-Clifford Generators and Operators
A slight generalization of the 1-D example we have discussed so far corresponds to a finite, discrete distribution (dv), that is, codomain of jumps  ↑ whose section  ↑ 1 is discrete and finite. A particular interesting case corresponds to this cone being generated by the canonical basis  (equation 50) of the generator algebra , that is, It corresponds indeed to a joint analysis or simulation of the operators e i forming this basis. The Figure 7 corresponds to this case for the pseudo-quaternion algebra. Equation (90) is easy to generalize to:

Earth and Space Science
10.1029/2019EA000608 where each half space (e i ) (equation 91) contains the vector Q v (e i )e i and is limited by the "hyperplane" e ⟂ i orthogonal to e i . Their intersection defines the (hyper) cube () whose edges are the semiaxes Q v (e i )e i , still due to equation (91). Tchiguirinskaia (2015, 2018) suggested such a (hyper) cube, but we presently derived it in a rigorous manner taking fully into account the signature of the quadratic form Q v . For instance, each of the four quarts of the plane generated by a pair of vectors (e i , e ) correspond to a configuration ({e i , e }) for a given signature(Q v (e i ), Q v (e )) and the 3-D figure 8 of Schertzer and Tchiguirinskaia (2018) presents ({e i , e , e k }) for the signature (+1, +1, +1).
Due to the orthogonality of the vector basis , equation (87), a cascade based on it reduces equation (87) to a series (e i ∈ ) of scalar equations all identical to equation (77): with the exception of the prefactor (e i ) that remains related to C 1 (e i ) through equation (81). However, it does not mean that these equations are statistically independent.
Having clarified the classical Lévy stable vector case, we can address the more general and more involved case that corresponds to and d and d(n) being matrices (Schertzer et al., 1999(Schertzer et al., , 2001. It was shown on Hilbert spaces that there is only a spectral constraint on the operator and that is which is a straightforward extension of the scalar case constraint (0 < ≤ 2). This is immediate for a diagonalizable ; otherwise, the deficiency of some generalized eigenspaces of −1 only introduces logarithmic corrections, as seen with the help of the Jordan form of −1 . Therefore, they do not change the spectral constraint based on power laws.
The fact that the quadratic form Q v of  is nondefinite does not either alter this spectral constraint. It is therefore surprising that the diagonalizable case has been so far disregarded, whereas it enables in a straightforward manner to have various stability indices, one for each eigenspace of , whereas the classical definition of stable vectors was limited to a unique stability index. For this case, it seems that all along the chain of derivations from the Poisson sum of jumps Y 's (equation 77) to the generator characteristic function (equation 96), one only needs to split the space in the direct sum of eigenspaces E i and to replace the scalar by the corresponding eigenvalue i . For instance, this would generalize equation (87) as follows: This generalization has to be further investigated.

Conclusions and Prospects
As mentioned in the abstract, this paper was motivated by a relative dissatisfaction of the present status of multifractals after three decades of their existence and a century for cascades. It is indeed ironic that we are only beginning to address multifractal vectors after such a long time devoted to scalar-valued multifractals, whereas Richardson's quatrain addressed the wind velocity. We have therefore argued that it is more than timely to bridge the gap between the potentials of multifractals and their actual application. After recalling some basic properties of the scalar framework, mostly with the help of the simplistic but pedagogic example of discrete cascades (section 3), we considered the responses of a field and a (mathematical) measure or generalized function to a scale change T of scale ratio , which are, respectively, given by the pullback transform T * and push forward T * , of T (section 4). We emphasized the fact that the group property of T just propagates to both transforms, and we gave examples of them in relation with scalar scaling. We also 10.1029/2019EA000608 recalled that a scale change T can be strongly isotropic (section 5). Considering noncommutative symmetry generators, we pointed out the interest of Clifford algebra as algebra of cascade generators (section 6). In section 7, the definition of directional statistical moments enabled us to generalize the definition of scalar multiscaling to operators: the (Laplace) second characteristic function of the generator needs only to have a logarithmic divergence with the scale ratio . We first gave the example of Gauss-Clifford operators, whose generators are Gaussian vectors of a Clifford algebra. The logarithmic divergence can be obtained with the help of a fractional integration and yields the moment scaling function of the Gauss-Clifford operators (equation (68)). Recalling that Gaussian laws are a very particular case of Lévy stable laws, whose main parameter is the stability index ∈ (0, 2), we considered the case of "strict" Lévy-Clifford operators ( < 2) that present the difficulty to be generated by jumps. This difficulty is overcome with the introduction of the codomain of jumps  ↑ and the domain of finite statistic orders  > because contrary to the 1-D case, they are not just the opposite of each other. The recognition of their structure of cones is particularly helpful, and both of them are indispensable to derive and concisely present the moment scaling function of the Lévy-Clifford operators (equations 96 and 101).
Foreseen developments are threefold. At first, the present results paved the way to physically based multifractal multivariate analysis. At the same time, they also require to extend the classical multifractal analysis techniques to this new framework, for example, how to extend to multifractal operators the Double Trace Moment technique that has been widely used for scalar multifractal fields. In particular, if the stability index of Lévy-Clifford generator, which is also the multifractality index of the resulting Lévy-Clifford operator, is a matrix instead of a scalar. Second, empirical analysis of multifractal operators must be actively developed to give some guidance when theoretical choices are difficult or even impossible. Indeed, the number of parameters can drastically increase with multifractal operators. Third, but not the least, it is indispensable to pursue the formalism developments to bring multifractal operators to operational level of simulations, as well as to consider extensions to processes that do not fall in the basin of application of the Lévy-Clifford algebra, despite their universal statistical and robust algebraic properties. In this direction, the relation between the Kiling form and the Clifford quadratic form should be assessed.