Laurent Baratchart
email: laurent.baratchart@inria.fr

Pei Dang
email: pdang@must.edu.mo.

Tao Qian

Hardy-Hodge Decomposition of Vector Fields in R n

We prove that a IR n+1 -valued vector field on IR n is the sum of the traces of two harmonic gradients, one in each component of IR n+1 \ IR n , and of a IR n -valued divergence free vector field. We apply this to the description of vanishing potentials in divergence form. The results are stated in terms of Clifford Hardy spaces, the structure of which is important for our study.

Introduction

Decomposing a complex function on the line as the sum of the traces of two holomorphic functions, one in each half plane cut out by the line, is a classical topic from complex analysis that lies at the root of many developments in harmonic analysis. Indeed, such a decomposition features the Hilbert transformation which is the prototype of a Calderòn-Zygmund operator, whose C 1,α and L p boundedness was historically the starting point of elliptic regularity theory [START_REF] Garnett | Bounded Analyic Functions[END_REF][START_REF] Muskhelishvili | Singular Integral Equations: Boundary Problems of Function Theory and Their Application to Mathematical Physics[END_REF][START_REF] Stein | Singular integrals and differetiability properties of functions[END_REF][START_REF] Torchinsky | Real variable methods in Harmonic Analysis[END_REF]. This decomposition is also a cornerstone of solutions to Riemann-Hilbert problems, which are especially meaningful in spectral theory [START_REF] Fokas | A unified transform method for solving linear and certain nonlinear PDEs[END_REF] and provided in recent years striking advances in the theory of orthogonal polynomials [START_REF] Deift | Orthogonal Polynomials and Random Matrices: a Riemann-Hilbert Approach[END_REF]. Moreover, it is instrumental for defining and studying Hankel and Toeplitz operators, which play a fundamental role in complex approximation and were successfully applied to issues of basic importance in control and signal analysis [START_REF] Nikolskii | Operators, Functions, and Systems: An Easy Reading, Part I&II[END_REF][START_REF] Peller | Hankel Operators and their Applications[END_REF][START_REF] Parfenov | Estimates of the singular numbers of a Carleson operator[END_REF][START_REF] Cohen | Time-Frequency Analysis[END_REF]. Besides, in a Hilbertain context, the decomposition was used to obtain sparse representations of analytic signals of scalar-valued signals in various classical contexts ( [START_REF] Qian | Adaptive Fourier Series-A Variation of Greedy Algorithm[END_REF] and subsequent papers by these authors).

Specifically, given a complex-valued function f ∈ L p (R), 1 < p < ∞, one has

f = f + + f -, (1.1)
where

f ± (x) = lim y→0± ±1 2πi ∞ -∞ f (t) t -(x + iy) dt,
where f ± are, respectively, non-tangential boundary limit functions of holomorphic functions of one complex variable in, respectively, the Hardy spaces H p (C ± ) of the upper and lower half planes. The Hardy space functions are given by

f ± (z) = ±1 2πi ∞ -∞ f (t)
t -z dt, z = x + iy, ±y > 0, see [START_REF] Garnett | Bounded Analyic Functions[END_REF]. Now, under the standard identification C ∼ IR 2 , a holomorphic function may be regarded as the gradient of a harmonic function, and this way (1.1) says that a IR 2 -valued vector field on IR of L p class may be decomposed as the sum of the traces of two harmonic gradients, arising from harmonic functions in the upper and lower half plane respectively. The question that we raise in this paper is whether such a decomposition is possible in higher dimension, namely whether a vector field in L p (IR n , IR n+1 ) is the sum of the traces of two gradients of functions harmonic in the two half-spaces cut out by IR n in IR n+1 . The answer is no in general, but the next best thing is that a decomposition becomes possible if a third summand is allowed, which takes the form of a divergence free vector field tangent to IR n . This fact was observed in [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF] when n = 2, and used to characterize silent magnetization distributions on a plane. We presently carry this decomposition over to every n. When projected onto IR n , it yields back the classical Hodge decomposition of a L p tangent vector field on IR n as the sum of a gradient and of a divergence free component. This why we call our decomposition of L p (IR n , IR n+1 ) vector fields the Hardy-Hodge decomposition.

Formally the decomposition can be surmised from Hodge theory for 1-currents supported on a hypersurface in ambient space [START_REF] Harris | Iterated integrals and cycles on algebraic manifolds[END_REF]Sec. 2.8], but the estimates needed to control L pnorms of the objects involved pertain to the Calderòn-Zygmund theory. In this connection, it would be pedantic to introduce currents to speak of vector fields on linear submanifolds, but it is convenient to use the formalism of Clifford analysis, which provides us with a substitute for complex variables and is well adapted to handle higher dimensional singular integrals. In fact, Clifford analysis is also suited to extend the result to vector fields on more general sub manifolds, although such a generalization lies beyond the scope of the present paper.

The latter is organized as follows. In Section 2 we recall some basic facts from Clifford analysis and Clifford Hardy spaces, most of which can be found in [START_REF] Gilbert | Clifford Algebras and Dirac Operators in Harmonic Analysis[END_REF], and we study the structure of boundary function in detail, along with density properties of rational-like functions. In Sections 3 and 4, we prove the Hardy-Hodge decomposition and some variants thereof. Finally, in Section 5, we discuss an application to non-uniqueness for inverse potential problems in divergence form.

Preliminaries

Let n ≥ 3 be an integer and Φ be either the real field IR or the complex field C. Hereafter, we put L p (IR n , E) for the familiar Lebesgue space of functions on IR n with values in a Banach space E (typically E = Φ m ) whose norm to the p-th power is integrable, and we often write L p (IR n ) for simplicity if E is understood from the context.

We adopt standard notations in Clifford analysis, see [START_REF] Gilbert | Clifford Algebras and Dirac Operators in Harmonic Analysis[END_REF]. In particular, we put Cl(n, Φ) to denote the Clifford algebra generated over Φ by e 1 , ..., e n with e 0 = 1 and e i e j + e j e i = -2δ ij , i, j = 1, ..., n. We indicate with S the collection of subsets of {1, • • • , n}. Then, the elements of the canonical basis of Cl(n, Φ), viewed as a vector space over Φ, are denoted as e S , S ∈ S, where e S = e j

1 • • • e j k if S = {1 ≤ j 1 < • • • < j k ≤ n}.
A generic member of Cl(n, Φ) can thus be written as x = Σ S∈S x S e S with x S ∈ Φ. When S is empty, we write e ∅ = e 0 = 1. The conjugate of x, denoted as x, is defined to be Σ S∈S (-1) |S| x S e S , where |S| indicates the cardinality of S. By convention, a 0-form is a scalar. A k-form is a sum Σ S∈S k x S e S where S k indicates those members of S with cardinality k. Clearly, Cl(n, Φ) is a 2 n -dimensional linear space over Φ. A 1-form is also called a vector, denoted with underscore:

x = x 1 e 1 + • • • + x n e n .
Clifford vectors are identified with Euclidean vectors in IR n . The sum of a 0-form and a 1 form is called a para-vector, and if x is a paravector we let x be its vector part: x = x 0 e 0 + x. This is consistent with our previous notation for vectors. The norm of x ∈ Cl(n, Φ) is defined to be |x| = ( S∈S |x S | 2 ) 1/2 , which derives from the inner product x, y = S∈S x S y S . If both x, y are para-vectors, then their Clifford product xy = -x, y + x ∧ y, where the exterior product x ∧ y is a 2-form similar to the exterior product of differential forms from geometry:

x ∧ y = Σ j<k (x j y k -x k y j )e j e k , x = Σ S∈S x S e S , y = Σ S∈S y S e S . Let R n 1 = {x = x 0 + x | x 0 ∈ R, x ∈ R n }, R n 1,± = {x = x 0 + x | ± x 0 > 0, x ∈ R n }, and R n+1 ± = {x = x + x n+1 e n+1 | ± x n+1 > 0}, noting that R n 1 , R n 1,± are included in Cl(n, IR) while R n+1 ± is contained in IR n+1 ⊂ Cl(n + 1, IR) (via
the identification of Clifford vectors with Euclidean vectors). For k = 0, ..., n + 1, we introduce the partial derivatives ∂ k = ∂/∂ x k and subsequently we define

D 0 = ∂ 0 , D n = e 1 ∂ 1 + • • • e n ∂ n , D = D 0 + D n , D n+1 = e 1 ∂ 1 + • • • e n ∂ n + e n+1 ∂ n+1 . (2.2) A Cl(n, Φ) or Cl(n + 1, Φ)-valued function f such that Df = 0 (f D = 0) or D n+1 f = 0 (f D n+1 = 0) on an open set of Ω ⊂ IR n+1 is called left-monogenic (right-monogenic) on Ω.
By convention, coordinates in the case of D are denoted by x 0 , • • • , x n , whereas in the case of D n+1 they are written x 1 , • • • , x n+1 . If a function is both left-and right-monogenic, we call it two-sided-monogenic. Let us stress that, when applying the differential operators (2.2), the partials ∂ j commute with the e k but the e j do not, so that it generally matters whether the operator gets applied from the left or the right. Note that

(D 0 -D n )D = ∆ (resp. D 2 n+1 = -∆) where ∆ = Σ n j=0 ∂ 2 x j (resp. ∆ = Σ n j=1 ∂ 2 x j
) is the ordinary Laplacian. Therefore left or right monogenic functions have harmonic components, in particular they are real analytic on Ω and there is no difference being monogenic in the distributional or in the strong sense. When f is para-vector valued in Cl(n, Φ) and we write

f = f 0 + f 1 e 1 + • • • + f n e n ,
it is readily checked that Df = 0 if and only if

∂ 0 f 0 = Σ n j=1 ∂ j f j with ∂ j f k = ∂ k f j and ∂ 0 f j = -∂ j f 0 for 1 ≤ j < k ≤ n, (2.3)
and similarly when f D = 0. In particular, a vector or para-vector valued function which is left-monogenic must be also right-monogenic. In the same manner, when f is vector valued in

Cl(n + 1, Φ) with f = f 1 e 1 + • • • + f n+1 e n+1 , we have that D n+1 f = 0 if and only if Σ n+1 j=1 ∂ j f j = 0 and ∂ j f k -∂ k f j = 0 for 1 ≤ j < k ≤ n + 1, (2.4)
and the same if f D n+1 = 0. Thus, vector valued left-monogenic functions are right-monogenic. By (2.3), a para-vector valued function

f = f 0 + f 1 e 1 + • • • + f n e n ,
where the f j are real valued, is monogenic if and only if (-f 0 , f 1 , • • • , f n ) is a harmonic gradient, meaning that it is the gradient of a harmonic function. The components of a harmonic gradient are sometimes referred to as a conjugate harmonic system, or a Riesz system of functions, cf. [START_REF] Stein | Intrucuction to Fourier analysis on Euclidean spaces[END_REF]. When IR n 1 gets identified with IR n+1 , the fact that (-

f 0 , f 1 , • • • , f n ) is a harmonic gradient amounts to saying that f = f 0 -f 1 e 1 -• • • -f n e n is a harmonic gradient. Likewise, it follows from (2.4) that a vector valued function in Cl(n + 1, IR), say f = f 1 e 1 + • • • + f n+1 e n+1 is monogenic if and only if (f 1 , • • • , f n+1
) is a harmonic gradient. Identifying vectors in Cl(n + 1, IR) with IR n+1 , we simply say in this case that f is a harmonic gradient.

Let g : R n 1,± → R n 1 . For 1 < p < ∞, we say g belongs to the Hardy space

H p (R n 1,± , R n 1 ) if Dg = 0 in R n 1,± and
g p H p ± sup ±x 0 >0 R n |g(x 0 + x)| p dx < ∞. (2.5)
We refer to the above Hardy spaces as being of para-vector type, or also of inhomogeneous type. Thus, g ∈ H p (R n 1,± , R n 1 ) if and only if g is a harmonic gradient which moreover satisfies the p-norm boundedness (2.5) in the relevant half space. Equivalently, since each of the functions composing a conjugate harmonic system is harmonic, it follows from (2.5) and standard estimates on harmonic functions (see e.g. [

g H p ± ≤ M α g L p (IR n ) ≤ C α g H p ± .
Here, to each α > 0 and x ∈ IR n , the notation Γ α (x) stands for the cone

Γ α (x) = {y 0 + y ∈ IR n 1,± , |y -x| < α|y 0 |},
and the precise value of α is irrelevant except that the constants will depend on it. Likewise, for g : R n+1 ± → R n+1 ⊂ Cl(n + 1, IR) and 1 < p < ∞, we say that g belongs to the Hardy space

H p (R n+1 ± , R n+1 ) if D n+1 g = 0 in the half space IR n+1 ± and g H p ±,h = sup ±x n+1 >0 R n |g(x + x n+1 e n+1 )| p dx < ∞. (2.7)
Note that, in the previous equation, x refers to a vector in Cl(n, IR) viewed as a vector in Cl(n + 1, IR) whose (n + 1)-st component is zero. We refer to the latter Hardy spaces as being of vector type or homogeneous type, which is the reason for the subscript "h" in the notation for the norm. Thus, we have that

f = f 1 e 1 + • • • + f n+1 e n+1 lies in H p (R n+1 ± , R n+1 ) if and only if it is a harmonic gradient in R n+1 ±
satisfying the p-boundedness condition (2.7). The latter is again equivalent to the L p boundedness on R n of the nontangential maximal function

M α,h g(x)
sup

x+x n+1 e n+1 ∈Γ α,h (x) |g(x + x n+1 e n+1 )| (2.8)
where this time Γ α,h (x) = {y + y n+1 e n+1 ∈ IR n+1 ± , |y -x| < α|y n+1 |}. In fact the passage from non homogeneous to homogeneous Hardy spaces is rather mechanical, trading x 0 and 1 for x n+1 and e n+1 while changing g 0 into -g n+1 .

Next, recall the local Fatou theorem asserting that a harmonic function in IR n+1 + which is non-tangentially bounded at almost every point of a set G ⊂ IR n has a non-tangential limit at almost every point of G; here, nontangential refers to the fact that bounds and limits are seeked in cones Γ α,h (x) for arbitrary but fixed α > 0, see [START_REF] Stein | Intrucuction to Fourier analysis on Euclidean spaces[END_REF]Thm. 3.19]. In view of (2.5) and (2.7), it follows from the local Fatou theorem that each component of a Hardy function (resp. homogeneous Hardy function) has a nontangential limit at almost every point of IR n ⊂ Cl(n, IR) (resp. IR n × {0} ⊂ IR n+1 ⊂ Cl(n + 1, IR)). This defines boundary values for such functions. Now, it is an important and peculiar property of left or right monogenic functions that they can be recovered as Cauchy integrals of their boundary values, see [START_REF] Gilbert | Clifford Algebras and Dirac Operators in Harmonic Analysis[END_REF]Cor. 3.20] and [START_REF] Kou | The Paley-Wiener theorem in Rn with the Clifford analysis setting[END_REF][START_REF] Li | Clifford algebras, Fourier transforms, and singular convolution operators on Lipschitz surfaces[END_REF]. Specifically, let us discuss the case of left monogenic Cl(n, IR)-valued functions on IR n 1,+ as a prototypical example. If g is such a function and if M α g ∈ L p (IR n ) for some p ∈ (1, ∞), then g has nontangential limit a.e. on IR n since its components are harmonic functions with L p nontangential maximal function, and this nontangential limit clearly lies in L p (IR n ) because it is dominated by the nontangential maximal function. Then, denoting the nontangential limit by g again and letting ω n designate the volume of the unit sphere of dimension n,, we get that g = C + g, where

C + g(x) = 1 ω n R n y -x |y -x| n+1 (-e 0 )g(y)dy (2.9) = 1 ω n R n x 0 |x -y| n+1 g(y)dy + 1 ω n R n x 1 -y 1 |x -y| n+1 (-e 1 )g(y)dy + • • • + + 1 ω n R n x n -y n |x -y| n+1 (-e n )g(y)dy, x = x 0 e 0 + x 1 e 1 + • • • + x n e n ∈ R n 1,+ ,
with ω n the surface measure of the n-dimensional unit sphere. Here, the presence of -e 0 in the definition of C + is because -e 0 is the exterior unit normal to IR n 1,+ , see [START_REF] Gilbert | Clifford Algebras and Dirac Operators in Harmonic Analysis[END_REF][START_REF] Li | Clifford algebras, Fourier transforms, and singular convolution operators on Lipschitz surfaces[END_REF]. Conversely, if g is any Cl(n, IR)-valued function in L p (IR n ), the formula for C + g makes good sense and defines a left monogenic function in IR n 1,+ since the kernel is left monogenic (as is easily verified). Invoking the Plemelj formula [START_REF] Li | Clifford algebras, Fourier transforms, and singular convolution operators on Lipschitz surfaces[END_REF], we get at almost all points x ∈ IR n that there exist non-tangential limits of C + g denoted with a curly C + g, which are given by

C + g(x) = 1 2 [g(x) + Hg(x)] = 1 2 (I + H)g(x), (2.10)
where We note that our definition of H is consistent with the classical definition of the Hilbert transform in dimension 1, given e.g. in [START_REF] Bell | The Cauchy Transform, Potential theory and Conformal Mappings[END_REF]. For if we designate the latter with H 1 , we get when n = 1, with e 1 = -i, that

Hg(x) = 2 ω n lim ε→0 + |x-y|>ε x 1 -y 1 |x -y| n+1 (-e 1 )g(y)dy + • • • + 2 ω n lim ε→0 + |x-y|>ε x n -y n |x -y| n+1 (-e n )g(y)dy n k=1 (-e k )R k (g)(x), where R k (g)(x) = 2 ω n lim ε→0 + |x-y|>ε x k -y k |x -y| n+1 g(y)
Hg = -e 1 1 π p.v. 1 (•) * g = i 1 π p.v. 1 (•) * g iH 1 g,
where "p.v." indicates the principal value. Thus, in view of the well known relation H 2 1 = -I, we have that H 2 = (iH 1 ) 2 = I, which is consistent with the n-dimensional case. Similar considerations apply to functions in H p (R n 1,-, R n 1 ), only trading C + for its opposite (because the outer normal to IR n 1,-is e 0 and not -e 0 ), and letting x now range over R n 1,-. Letting this time C -g stand for the nontangential limit on IR n , this results in the Plemelj formula:

C -g(x) = 1 2 (I -H)g(x).
(2.11)

In the sequel we denote by Sc{x} the scalar part of a Clifford number x ∈ Cl(n, Φ), which is the 0-form of x, and by Nsc{x} the non-scalar part of x, which is the sum of all the k-forms of x, k = 1, ..., 2 n .

We pointed out already that each function in

H p (R n 1,± , R n 1
) is naturally associated to the conjugate of a harmonic gradient. In the next lemma, we identify the latter as being the gradient of a Newton potential, and we describe the boundary values of Hardy functions.

Lemma 2.1. If g = g 0 e 0 + • • • + g n e n ∈ H p (R n 1,± , R n 1 )
, where each g k is real valued, then its non-tangential limit on IR n , still denoted as g, satisfies g = (I ± H)g 0 , that is,

g k = ∓R k g 0 for k = 1, ..., n. Conversely, each function on IR n of the form (I ± H)ϕ, with ϕ ∈ L p (IR n , IR), is the nontangential limit of a function in H p (R n 1,± , R n 1 )
and the Hardy norm is equivalent to

ϕ L p (IR n ) . Moreover, the conjugate of each g ∈ H p (R n 1,± , IR n 1 )
is the harmonic gradient of the Newton potential of ∓2/(n -1) times its scalar part g 0 , namely:

g(x) = ∂ 0 N g 0 (x)-Σ n k=1 ∂ k N g 0 (x)e k , N g 0 (x) ∓ 2 (n -1)ω n IR n g 0 (y) |x -y| n-1 dy, x ∈ IR n 1,± . (2.12)
Proof. By the Cauchy-Clifford formula, the function g in the upper-half space may be expressed as the Cauchy integral over R n of its boundary limit function. Therefore the Plemelj theorem implies that (1/2)(I ± H)g = g, hence ±Hg = g. By comparing the terms of degree 1 on the left-and the right-hands of the above identity we get ∓e k R k g 0 = e k g k for k = 1, ..., n, ensuing that g = (I ± H)g 0 . Invoking the Plemelj theorem again, this implies that the Cauchy integral of 2g 0 is a left monogenic function whose nontangential boundary value is g. Hence, by the Cauchy formula, this function must coincide with g on IR n 1,± . Now, when identifying IR n 1 with IR n+1 , the conjugate of the Cauchy kernel is -1/(n -1) times the gradient of the Newton kernel 1/|x -.| n-1 , and taking gradient commutes with the integration, as the integrand does not have singularity, so we conclude that the conjugate of g is the gradient of the Newton potential of ∓2g 0 /(n -1) (the ∓ arises because of the presence of ∓e 0 in the Cauchy integral). Conversely, by (2.10), a function on IR n of the form (I ± H)ϕ, where ϕ is scalar valued in L p (IR n ), is the nontangential limit of the Cauchy integral C ± (2ϕ) which is indeed para-vector valued an monogenic. To see that it lies in H p (R n 1,± , IR n 1 ), observe from what precedes that its conjugate is the gradient of the Newton potential N ϕ of ∓2ϕ/(n -1). In particular, by inspection of formula (2.12) (where g 0 is set to ϕ), we find that Sc{(C ± ϕ)(x 0 , x)} = (P x 0 * ϕ)(x) at every (x 0 , x) ∈ IR n 1,± , where the symbol " * " indicates convolution and P x 0 is the Poisson kernel at level x 0 :

P x 0 (x) = 2 ω n |x 0 | (x 2 0 + |x| 2 ) (n+1)/2 , x ∈ IR n . (2.13)
Since P x 0 has unit norm in L 1 (IR n ) for all x 0 , it follows that

Sc (C ± ϕ)(x 0 , .) L p (IR n ) ≤ ϕ L p (IR n ) , ∀x 0 , (2.14)
implying that the scalar part of C ± ϕ meets the p-boundedness condition (2.5). To show that the vector part also satisfies this condition, let us work on IR n 1,+ as the argument on IR n 1,-is similar. Fix z 0 > 0 and consider the para-vector valued function on IR n 1,+ given by F (x 0 , x) = (C + ϕ)(z 0 +x 0 , x). Clearly it is monogenic, and we get upon applying Hölder's inequality to (2.9) that F H p + ≤ c for some constant c = c(z 0 ), see definition (2.5). Hence F ∈ H p (R n 1,± , IR n 1 ), and it is obvious that its nontangential limit on IR n is C + ϕ(z 0 , .). thus, by the previous part of the proof, it holds that

(C + ϕ)(x 0 , .) = (I + H)Sc{F (0, .)} = (I + H)Sc (C + ϕ)(x 0 , .) .
In view of the definition of H (cf. equation following (2.10)) and the L p -boundedness of Riesz transforms, we now deduce from (2.14) that C + ϕ satisfies (2.5), as desired. We also proved that

C + ϕ H p + ≤ c ϕ L p (IR n )
for some constant c independent of ϕ, and since ϕ L p (IR n ) is obviously less than the L p norm of the non-tangential maximal function which itself is equivalent to the Hardy norm, as pointed out after (2.6), the proof is complete.

When dealing with functions in H p (R n+1 ± , R n+1 ), the Cauchy formula is the same except that ∓e 0 gets replaced by ∓e n (the outer normal to IR n+1 ± ), and in the Plemelj formula H is changed into He n+1 = Σ n k=1 (-e k e n+1 )R k , see [START_REF] Gilbert | Clifford Algebras and Dirac Operators in Harmonic Analysis[END_REF]Ch. 2,Sec. 5]. An argument analogous to the previous one shows that

f = f 1 e 1 + • • • + f n+1 e n+1 lies in H p (R n+1 ± , R n+1 ) if and only if f k = ±R k f n+1 for k = 1, ..., n. and that f is D n+1 N ,
where N is the Newton potential of ∓2g n+1 /(n -1).

Observe that Lemma 2.1 and its analog for homogeneous spaces entail that

H p (IR n 1,± , IR n 1 ) (resp. H p (IR n+1 ± , IR n+1
)) is a Banach space isomorphic to L p (IR n , IR), with norm equivalent to the L p -norm of the trace of the scalar part (resp. n + 1-st component). Observe also from this, since L q (IR n , IR)

∩ L p (IR n , IR) is dense in L p (IR n , IR) for any q ∈ (1, ∞), that H p (R n+1 ± , R n+1 ) ∩ H q (R n+1 ± , R n+1 ) is dense in H p (R n+1 ± , R n+1
). Lemma 2.1 easily implies a result which is of interest in its own right and parallels the density of rational functions in holomorphic Hardy spaces of index p ∈ (1, ∞) on the half-plane [START_REF] Garnett | Bounded Analyic Functions[END_REF]. Note that rational functions with simple poles are conjugate of gradients of logarithmic potentials of discrete measures with finite support. In the present context, analogs of rational functions with simple poles are conjugates of gradients of Newton potentials of discrete measures with finite support. Specifically, if we let

R x (y) = ∇ y 1 ω n |x -y| n-1 = n -1 ω n x 0 -y 0 |x -y| n+1 - n j=1 x j -y j |x -y| n+1 e j , y ∈ IR n 1,+ , x ∈ IR n 1,-, then R x ∈ H p (IR n 1,+ , IR n 1 )
L p (IR n , IR) h → C + h -Σ n k=1 (R k h)e k ∈ H p (IR n 1,+ , IR n 1 ). (2.15)
The inverse image of R x under this isomorphism is (1-n)/2 times the Poisson kernel P x 0 (x-y) defined in (2.13). Thus, by the Hahn-Banach theorem, the asserted density is equivalent to the fact that no nonzero function in L p (IR n ), IR), with 1/p + 1/p = 1, can have vanishing Poisson integral. This, however, drops out immediately from the property that the Poisson kernel is an approximate identity, thereby achieving the proof. When saying that a vector-valued function

f = f 1 e 1 + • • • + f n e n on R n is divergence free, we mean that divf = n k=1 ∂ k f k = 0
. This is to be understood in the generalized function sense that amounts to the relation n k=1 R k f k = 0 or, equivalently, n k=1 ξ k fk (ξ) = 0, through the inverse Fourier transformation, again to be understood in the generalized function sense if p > 2 so that the Fourier transform is really a distribution. The space of vector valued divergence free maps in L p (IR n , IR n ) is a closed subset thereof and thus a Banach space in its own right that we denote by D p (IR n ). Though initially defined on IR n only, a divergence free vector field extends naturally to IR n 1,± (resp. IR n+1 ± ) into a IR n valued map F which is independent of x 0 (resp. x n+1 ). This function needs not be monogenic, but it satisfies Sc DF = 0 (resp. Sc D n+1 F = 0).

Hardy-Hodge Decomposition of Para-Vector-Valued functions in

L p (R n ) Theorem 3.1.
Let f be a para-vector valued function in L p (R n , Cl(n, Φ)), 1 < p < ∞. Then f is uniquely decomposed as f = f + + f -+ f 0 , all in L p (R n ), such that f ± are para-vector-valued, being the non-tangential boundary limits of some two functions in, respectively, H p (R n 1,± , R n 1 ), and f 0 is vector-valued and divergence free. Moreover, for all p in the indicated range the decomposition is unique and, for p = 2,

f 2 = f + 2 + f -2 + f 0 2 . (3.16)
In fact, this decomposition induces a topological direct sum:

L p (IR n , IR n 1 ) = H p (R n 1,+ , R n 1 ) ⊕ H p (R n 1,-, R n 1 ) ⊕ D p (IR n ). (3.17) Proof Let f (x) = n k=0 f k (x)e k be in L p (R n ), where f 0 (x), f 1 (x), • • • , f n (x) are scalar-valued and x = x 1 e 1 + • • • + x n e n ∈ R n .
By (2.10) and (2.11), since H 2 = I, we have that C ± is a projection:

1 2 (I + H)f (x) = [ 1 2 (I + H)] 2 f (x), 1 2 (I -H)f (x) = [ 1 2 (I -H)] 2 f (x).
Then we have

f (x) = 1 2 (I + H)f (x) + 1 2 (I -H)f (x) = [ 1 2 (I + H)] 2 f (x) + [ 1 2 (I -H)] 2 f (x) = 1 2 (I + H) Sc{ 1 2 (I + H)f (x)} + Nsc{ 1 2 (I + H)f (x)} + 1 2 (I -H) Sc{ 1 2 (I -H)f (x)} + Nsc{ 1 2 (I -H)f (x)} = 1 2 (I + H)[Sc{ 1 2 (I + H)f (x)}] + 1 2 (I -H)[Sc{ 1 2 (I -H)f (x)}] + 1 2 (I + H)[Nsc{ 1 2 (I + H)f (x)}] + 1 2 (I -H)[Nsc{ 1 2 (I -H)f (x)}]. (3.18)
Consider the function given by the last line of the above chain of equalities, viz.

1 2 (I + H)[Nsc{ 1 2 (I + H)f (x)}] + 1 2 (I -H)[Nsc{ 1 2 (I -H)f (x)}]. (3.19)
It can be computed directly through

1 2 (I + H)[Nsc{ 1 2 (I + H)f (x)}] = 1 4 (I + H)[ n k=1 f k e k + Nsc{Hf }] = 1 4 n k=1 f k e k + Nsc{Hf } + n k=1 H[f k e k ] + H[Nsc{Hf }] ,
and

1 2 (I -H)[Nsc{ 1 2 (I -H)f (x)}] = 1 4 (I -H)[ n k=1 f k e k -Nsc{Hf }] = 1 4 n k=1 f k e k -Nsc{Hf } - n k=1 H[f k e k ]] + H[Nsc{Hf }] .
By adding these relations together, we have that

1 2 (I + H)[Nsc{ 1 2 (I + H)f (x)}] + 1 2 (I -H)[Nsc{ 1 2 (I -H)f (x)}] = 1 2 n k=1 f k e k + H[Nsc{Hf }] . (3.20)
As f is para-vector valued, relation (3.18) and the fact that H maps scalar valued functins to vector valued functions together imply that the quantity in (3.20) is a paravector, therefore H[Nsc{Hf }] is a para-vector. Now we work out its expression. Since

Hf = n k=1 R k (f k ) - n k=1 R k (f 0 )e k + H ∧ f , where f = n k=1 f k e k , we have that Nsc{Hf } = - n k=1 R k (f 0 )e k + H ∧ f .
Consequently, since we need only collect terms of degree 0 and 1, and because Riesz transforms commute, we obtain:

H[Nsc{Hf }] = (- n k=1 R 2 k )f 0 + n k=1 ( l =k -R 2 l )f k + R k l =k R l f l e k = f 0 + n k=1 ( l =k -R 2 l )f k + R k l =k R l f l e k .
Substituting back into (3.20), we get that

1 2 (I + H)[Nsc{ 1 2 (I + H)f (x)}] + 1 2 (I -H)[Nsc{ 1 2 (I -H)f (x)}] = 1 2 f (x) + n k=1 ( l =k -R 2 l )f k + R k l =k R l f l e k . .
Therefore, by (3.18),

f (x) = 1 2 (I + H) Sc{ 1 2 (I + H)f (x)} + 1 2 (I -H) Sc{ 1 2 (I -H)f (x)} + 1 2 f (x) + n k=1 ( l =k -R 2 l )f k + R k l =k R l f l e k .
Finally,

f (x) = (I + H) Sc{ 1 2 (I + H)f (x)} + (I -H) Sc{ 1 2 (I -H)f (x)} + + n k=1 ( l =k -R 2 l )f k + R k l =k R l f l e k . ( 3
( l =k -R 2 l )f k + R k l =k R l f l e k
is divergence-free. For this, by the last remark of the last section, it suffices to show that

n k=1 R k ( l =k -R 2 l )f k + R k l =k R l f l = 0. (3.22)
The above, however, is obvious since Riesz transformations commute, and thus we obtain the desired decomposition f (x) = f + + f -+ f 0 , where

f + = (I + H)[Sc{ 1 2 (I + H)f (x)}] ∈ H p (R n 1,+ , IR n 1 ), (3.23) f -= (I -H)[Sc{ 1 2 (I -H)f (x)}] ∈ H p (R n 1,-, IR n 1 ), (3.24)
and

f 0 = n k=1 ( l =k -R 2 l )f k + R k l =k R l f l e k (3.25)
is divergence free. Next we prove uniqueness. This is equivalent to showing that if we have a decomposition of the zero function: 0 = f + + f -+ f 0 , then it must be that f + = f -= f 0 = 0. Indeed, in that case we may write

0 = (I + H)(f + + f -+ f 0 ) = 2f + + (I + H)f 0 = 2f + + f 0 + Hf 0 ,
where we used that f + = Hf + by (2.10) since C + f + = f + , and also that Hf -= -f -by (2.11) since C -f -= f -. Note that since f 0 is divergence free, the scalar part of Hf 0 is zero, and thus only the 2-form part of Hf 0 is possibly non-zero. However, the last equality shows that the 2-form part also has to be zero, because all the other terms are para-vectors. We thus conclude that 2f + + f 0 = 0. The same reasoning gives us 2f -+ f 0 = 0. These together yield f + = f -. By applying I + H to both sides we get f + = 0 hence also f -= 0, and consequently f 0 = 0. This establishes uniqueness and shows that (3.17) holds as a direct sum. In addition, since f ± and f 0 are continuous functions of f in L p (IR n , Cl(n, Φ)) by (3.23), (3.24), (3.25) and the L p continuity of the Riesz transformations, we see that the projections in (3.17) are continuous, hence the sum is topological by the open mapping theorem. Finally, when p = 2, we show the Pythagora type relation

f 2 = f + 2 + f -2 + f 0 2 . (3.26)
First, since f is para-vector-valued, we obviously have that

f 2 = R 2 f f dx.
Hence, to prove (3.26), it suffices to establish the following orthogonality relations:

R n f + f -dx = R n f -f + dx = 0, (3.27)
and

R n (f + f 0 + f 0 f + )dx = R n (f -f 0 + f 0 f -)dx = 0. (3.28) Let us show (3.27). Recall from Lemma 2.1 that if g ∈ H 2 (R n 1,± , R n 1 )
, then g = (I ± H)g 0 and consequently, taking Fourier transforms (the Fourier transform of a vector valued function is computed componentwise), we get that ĝ = 2χ ± ĝ0 , where χ ± are multipliers for the Hardy space projections:

χ ± (ξ) = 1 2 (1 ± i ξ |ξ| )
, that satisfy χ 2 ± = χ ± and χ + + χ -= 1 as well as χ + χ -= χ -χ + = 0. Here we used the expression for the multiplier of R k in the Fourier domain, see [START_REF] Stein | Singular integrals and differetiability properties of functions[END_REF][START_REF] Li | Clifford algebras, Fourier transforms, and singular convolution operators on Lipschitz surfaces[END_REF]. Applying these remarks to g = f ± , and using Parseval's Theorem, we have that

R 2 f + f -dx = R 2 (f + ) ∧ (ξ)(f -) ∧ (ξ)dξ = 4 R 2 χ + (ξ)χ -(ξ) f + 0 f - 0 dξ = 0,
where we used the relation χ + (ξ)χ -(ξ) = 0 for all ξ. The proof of the second equality relation in (3.27) is similar. Now, let us show that equality (3.28) holds. Indeed,

R n (f + f 0 + f 0 f + )dx = 2Sc R n f + f 0 dx = Sc R n 1 + i ξ |ξ| f 0 (ξ) f 0 (ξ)dξ = Sc R n f 0 (ξ) iξ |ξ| f 0 (ξ)dξ = 0,
where the last equality used the relation (iξ/|ξ|) f 0 (ξ) = 0, the latter being a consequence of the fact that f 0 is divergence free. The proof is complete.

Remark If, alternatively, we use the scalar product

f, g = Sc R n f gdx, then indeed the decomposition f = f + + f -+ f 0 is orthogonal. Moreover, since we observed after the proof of Lemma 2.1 that H 2 (IR n 1,± , IR n 1 ) ∩ H p (IR n 1,± , IR n 1 ) is dense in H p (IR n 1,± , IR n 1 )
, we deduce that (3.27) holds as soon as

f + ∈ H p (IR n 1,+ , IR n 1 ) and f -∈ H p (IR n 1,-, IR n 1 ) with 1/p + 1/p = 1: . R n f + f -dx = R n f -f + dx = 0, f + ∈ H p (IR n 1,+ , IR n 1 ), f -∈ H p (IR n 1,-, IR n 1 ). (3.29)
Likewise, (3.28) generalizes to

R n (f + f 0 + f 0 f + )dx = R n (f -f 0 + f 0 f -)dx = 0, f ± ∈ H p (IR n 1,± , IR n 1 ), f 0 ∈ D p (IR n ). (3.30)
A few comments are in order:

• If (3.17) gets projected onto the last n components, and since (±R 1 h, • • • , ±R n h) is a gradient vector field on IR n (i.e. the gradient of the trace of a solution to the Neumann problem on IR n+1 ± with inner normal derivative h a.e. on IR n ) we recover the classical Helmoltz-Hodge decomposition of vector fields from L p (IR n , IR n ) into the sum of a rotational-free and a divergence-free vector field [START_REF] Iwaniec | Geometric Function Theory and Non-linear Analysis[END_REF].

• Decomposition (3.17) generalizes to higher dimensions the standard decomposition of a complex valued function in L p (IR) into the sum of a function belonging to the holomorphic Hardy space H p (IR 2 + ) and a function belonging to the holomorphic Hardy space H p (IR 2 -). The difference in dimension bigger than 1 is that a divergence free term must be added, for in this case not every vector field is a gradient. Note, since R k and the divergence operator preserve realness, that Theorem 3.1 carries over to Clifford valued maps and Clifford Hardy spaces with complex coefficients.

Variations

Next we consider the homogeneous case on R n . We regard IR n as being the subspace IR n ×{0} of IR n+1 . When considering Clifford 1-forms as Euclidean vectors, it means that 1-forms in Cl(n, IR) get identified with 1-forms in Cl(n + 1, IR) whose coefficient of e n+1 is zero. Note that Cl(n, IR) can be viewed as the subalgebra of Cl(n + 1, IR) generated by e 1 , • • • , e n . Theorem 4.1. Let f ∈ L p (R n , Cl(n + 1, IR)) be an (n + 1)-vector-valued function, 1 < p < ∞. Then f is uniquely decomposed as f = f + + f -+ f 0 , all in L p (R n ), such that f ± are the non-tangential boundary limits of some two functions in, respectively, H p (R n+1 ± , R n+1 ), while f 0 is vector-valued in Cl(n, IR) and divergence free. For all p in the indicated range, the decomposition is unique and topological. Moreover, for p = 2,

f 2 = f + 2 + f -2 + f 0 2 .
By factorizing out e n+1 one can, in particular, reduce the proof of Theorem 4.1 to that of Theorem 3.1. More precisely, noting that e -1 n+1 = -e n+1 , we use the relation

n+1 k=1 f k e k = [ n+1 k=1 f k e k e -1 n+1 ]e n+1 .
For the Cauchy kernels in the two settings, one has This corresponds to the relation between the two Dirac operators:

n+1 k=1 ∂ k e k = ( n+1 k=1
∂ k e k e -1 n+1 )e n+1 .

Letting ẽk = e k e -1 n+1 , k = 1, ..., n, one reduces the proof of Theorem 4.1 to that of Theorem 3.1 (compare [START_REF] Qian | Fourier analysis on starlike Lipschitz surfaces[END_REF]).

Theorem 4.1 can alternatively be rewritten without mentioning Clifford analysis as Theorem 4.2. Every vector field f ∈ L p (R n , R n+1 ), 1 < p < ∞, may be uniquely decomposed as f = f + + f -+ f 0 , where f ± ∈ L p (R n , R n+1 ) are, respectively, the non-tangential boundary limits of some harmonic gradients on R n+1 ± , which satisfy (2.7), while f 0 ∈ L p (IR n , IR n ) is divergence free. The decomposition is topoligical and for p = 2 there holds

f 2 = f + 2 + f -2 + f 0 2 .
From the analog of Lemma 2.1 for homogeneous Hardy spaces (see the discussion after the proof of that lemma), we know that boundary limits of harmonic gradients on R n+1 ± satisfying (2.7) are those members of L p (IR n , IR n+1 ) of the type (±R 1 h, • • • , ±R n h, h) with h ∈ L p (IR n , IR), where all the "+" signs and, respectively, all the "-" signs are taken. In such form, Theorem 4.2 was proven in [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF] when n = 2, also for more general function spaces.

It is worth contrasting Theorem 3.1 and Theorem 4.1 with their quaternionic counterparts. The space H of real quaternions consists of numbers q = q 0 +q, with q = q 1 e 1 +q 2 e 2 +q 3 e 3 where q j ∈ IR and e 1 , e 2 are as before, but, additionally, e 3 = e 1 e 2 . We identify the linear space consisting of all q with the space R 3 , and we put H ± for those quaternions with, respectively, ±q 0 > 0. We say that a quaternionic valued function f is left quaternionic if (D 0 +D 3 )f = 0, see (2.2), but this time the relation e 1 e 2 = e 3 is taken into account. The definitions of quaternionic Hardy spaces H p (H ± , H) as spaces of left quaternionic functions in H ± meeting the analogs of (2.5) now run parallel to those for inhomogeneous Hardy spaces. 3 , H) are non-tangential boundary limit functions of some two functions in, respectively, H p (H ± , H). Moreover, for p = 2, f + and f -are orthogonal:

Theorem 4.3. Let f ∈ L p (R 3 , H), 1 < p < ∞. Then f is uniquely decomposed as f = f + + f - such that f ± ∈ L p (R
f 2 = f + 2 + f -2 .
Proof. This is an immediate consequence of the corresponding Plemelj formula, for functions obtained from the Cauchy formula are all quaternionic-valued.

Remark: the reason why the divergence free term f 0 can be omitted in Theorem 4.3 is the closedness of multiplication in the quaternionic field. At the same time, the interpretation of f ± as traces of harmonic gradients is lost. To see the difference in the analysis of H p (R 3 1,± , R 3 1 ) and H p (H ± , H), recall that if f (x 0 +x) = f 0 (x 0 +x)e 0 +f 1 (x 0 +x)e 1 +f 2 (x 0 +x)e 2 +f 3 (x 0 +x)e 3 lies in H p (R 3 1,+ , R

1 ), then      ∂f 0 ∂x 0 = ∂f 1 ∂x 1 + ∂f 2 ∂x 2 + ∂f 3 ∂x 3 , ∂f 0 ∂x i = -∂f i ∂x 0 , i = 1, 2, 3, ∂f i ∂x j = ∂f j ∂x i , i = 0, j = 0, i = j, 3
while f (q 0 + q) = f 0 (q 0 + q)e 0 + f 1 (q 0 + q)e 1 + f 2 (q 0 + q)e 2 + f 3 (q 0 + q)e 3 ∈ H p (H + , H) will imply

         ∂f 0 ∂x 0 = ∂f 1 ∂x 1 + ∂f 2 ∂x 2 + ∂f 3 ∂x 3 , ∂f 1 ∂x 2 = ∂f 3 ∂x 0 + ∂f 0 ∂x 3 + ∂f 2 ∂x 1 ∂f 2 ∂x 3 = ∂f 1 ∂x 0 + ∂f 0 ∂x 1 + ∂f 3 ∂x 2 ∂f 3 ∂x 1 = ∂f 2 ∂x 0 + ∂f 0 ∂x 2 + ∂f 1 ∂x 3 . ( 4

.32)

The system of equations (4.31) implies the system of equations (4.32), but not conversely.

Uniqueness issues for potentials in divergence form

In [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF], the Hardy-Hodge decomposition was introduced when n = 2 to characterize silent magnetizations with support in IR 2 . This issue can be recast as that of describing vanishing potentials in divergence form. Recall that the Newton potential of a distribution ϕ on IR n+1 is the convolution of ϕ with 1/(ω n |x| n-1 ), wherever it exists. The potential is said to be in divergence form if ϕ can be taken to be the divergence of some IR n+1 valued distribution ψ:

P div ψ (x) = 1 ω n 1 |x -y| n-1 div ψ(y) = - 1 ω n ∇ y 1 |x -y| n-1 • ψ(y), (5.33)
where ∇ y indicates the gradient with respect to the variable y, where the dot indicates Euclidean scalar product. When n = 2, in the quasi-static approximation to Maxwell's equations, (5.33) formally expresses the magnetic potential of the magnetization ψ (cf. [10, Section 5.9.C]). Those x (if any) for which this expression makes good sense depend of course on ψ. We shall be concerned with the case where ψ is supported on a hyperplane P and has L p density there. Specifically, if we write

P = {x ∈ IR n+1 , x • u = a} (5.34) for some u ∈ IR n+1 and a ∈ IR, it means that ψ = f ⊗ δ 0 (x • u -a) where f = (f 0 , • • • , f n ) ∈ L p (P, IR n+1
) and δ 0 indicates the Dirac mass at 0 . Then, (5.33) becomes

P div ψ (x) = n -1 ω n P f (y) • (x -y) |x -y| n+1 dy, (5.35)
which is well defined for all x / ∈ P, more generally for all x not in the support of f . Let

H ± = {x ∈ IR n+1 , ±(x • u -a) > 0}
denote the two half spaces whose union is IR n+1 \ P, the complement of P. The question that we raise is:

for which f does it happen that P div ψ (x) = 0 for all x ∈ H ± ?

From the physical viewpoint, it amounts when n = 1 to describe those magnetizations with L p density supported on a plane which are silent from one side of that plane, meaning that they generate no magnetic field in the corresponnding half-space. These cannot be detected by measuring devices and account a good deal for the ill-posedness of inverse magnetization problems [START_REF] Parker | Geophysical Inverse Theory[END_REF]. The result below gives an answer to the question in terms of the Hardy-Hodge decomposition, thereby generalizing to higher dimension results from [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF] for n = 2. It will be convenient to define the Clifford Hardy spaces H p (H ± , IR n 1 ) consisting of paravector valued monogenic functions g in H ± meeting the condition: where dm indicates the differential of Lebesgue measure. Just like in the case of inhomogeneous Hardy spaces H p (IR n 1,± , IR n 1 ), condition (5.36) may be replaced by the L p (P, IR)-boundedness of the nontangential maximal function, computed this time over cones with vertex on P and axis parallel to u. Functions in H p (H ± , IR n 1 ) have nontangential limits in L p (P, IR n 1 ) of which they are the Cauchy Clifford integral, and they can be identified with their nontangential limit. In fact, if R is any orientation preserving affine isometry mapping P to {0} × IR n , we have that f belongs to H p (H ± , IR n 1 ) if and only if f • R belongs to H p (IR n 1,± , IR n 1 ). As an extra piece of notation, we put u ⊥ for the vector space orthogonal to u in IR n+1 , which is the linear hyperplane parallel to P. Theorem 5.1. Let P ⊂ IR n+1 be a hyperplane defined by (5.34) and ψ ∈ L p (P, IR n+1 ) with 1 < p < ∞. Then the potential P divψ vanishes on H ± if and only if ψ is the sum of a member of H p (H±, IR n 1 ) and of a divergence free vector field tangent to P in L p (P, u ⊥ ). The potential P divψ vanishes on IR n+1 \ P (that is, on both H + and H -) if and only if it is a divergence free function in L p (P, u ⊥ ).

Proof. Because the statement is invariant under orientation preserving affine isometries, we may assume that P = {0} × IR n so that H ± = IR n 1,± . Let us single out the components of ψ as ψ 0 , • • • , ψ n and identify ψ with the para-vector valued function ψ = ψ 0 e 0 + ψ 1 e 1 + • • • + ψ n e n . Set ψ = ψ + + ψ -+ ψ 0 for the Hardy-Hodge decomposition from Theorem 3.1. For x ∈ IR n 1,± , it is easily checked that y → x -y/|x -y| n+1 lies in H q (IR n 1,± , IR n 1 ). Thus, it follows from (5.35), (3.29) and (3.30) that P div(ψ ∓ +ψ 0 ) ≡ 0 on IR n 1,± . Therefore, the assumption that P divψ ≡ 0 on IR n 1,± reduces to P divψ ∓ ≡ 0 on IR n 1,± . Now, comparing (5.35) and (2.9), we find this is equivalent to Sc C ± ψ ∓ (x) = 0, x ∈ IR n 1,± . which amounts to Sc C ∓ ψ ∓ (ξ) = 0, ξ ∈ IR n 1,∓ .

(5.37) Since ψ ∓ ∈ H p (IR n 1,∓ , IR n 1 ), we have by the Cauchy Clifford formula that C ∓ ψ ∓ (ξ) = ψ ∓ (ξ), therefore (5.35) means that Sc{ψ ∓ } vanishes on IR n 1,∓ and so does its nontangential limit on IR n . But we know from Lemma 2.1 that the L p (IR n )-norm of the nontangential limit of the scalar part is an equivalent norm on H p (IR n 1,∓ , IR n 1 ), hence ψ ∓ = 0. This proves the first assertion of the theorem. To establish the second assertion, observe from what precedes that if P div ψ = 0 both in IR n 1,+ and IR n 1,-, then ψ ± = 0, and thus ψ = ψ 0 is vector-valued and divergence free. The proof is complete.

  dy is the k-th Riesz transformation of g, k = 1, ..., n. As an operator, R k has multiplier -iκ k /|κ| in the Fourier domain (κ = (κ 1 , • • • , κ n ) denoting the Fourier variable) and it maps L p (IR n , IR) into itself for 1 < p < ∞, see [20, Ch. III, Secs. 4.3-4.4]. The operator H defines the Hilbert transformation in the present context, and it satisfies the relation H 2 = I, where I is the identity.

n k=1 y k

  e k -( n+1 k=1 x k e k ) | n k=1 y k e k -( n+1 k=1 x k e k )| n+1 = -e n+1 n k=1 y k e k e -1 n+1 -( n k=1 x k e k e -1 n+1 + x n+1 ) | n k=1 y k e k e -1 n+1 -( n k=1 x k e k e -1 n+1 + x n+1 )| n+1 .

  sup ±b>0 P b |g| p dm < ∞, P b {x ∈ IR n+1 , ±(x • u -a) = b},(5.36)

Acknowledgements

This work was supported by Macao Science and Technology Development Fund, MSAR. Ref. 045/2015/A2,The work was supported by Multi-Year Research Grant (MYRG) MYRG116(Y1-L3)-FST13-QT. Macao Government FDCT 098/2012.


Richard Malgat

Benjamin Gilles

David I W Levin

Matthieu Nesme

Franc ¸ois Faure

Multifarious Hierarchies of Mechanical Models for Artist Assigned Levels-of-Detail

Keywords: CR Categories: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling-Physically based modeling I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation; Physically Based Animation, Deformable Solid, Multiscale Continuum Mechanics

a) Coarse framebased simulation (b) Fine FEM simulation (c) Mixing a coarse frame-based simulation to a local FEM patch with our method

Figure 1: Using a hook to catch a coarsely discretized frame-based octopus is challenging due to the lack of deformation at the contact point (1a). One could resolve this problem with a high-resolution set of finite elements (1b) but, at the expense of runtime performance and with the introduction of spurious secondary motion. Multifarious Hierarchies of Mechanical Models allow an artist to add arbitrarily located detail simulations to the underlying coarse deformation models in order to produce the desired result. (1c).

Introduction

Physics-based animation has become an important arrow in the quiver of the visual effects practitioner -routinely used to produce creatures and scenes that appear plausible while simultaneously defying belief. In truth these effects depend on artists carefully sculpting their simulation results as much as they do on the simulations themselves. It is not uncommon for an artist to produce a highly scripted animation and then use a physical simulator to augment it with secondary motion (skin wrinkling, muscles flexing, metal bending). And, like a painter adding the finishing strokes to a canvas, a visual effects artist has an innate understanding of the scale, location and type of motion that should be added to the final scene.

Adapting the level-of-detail of a simulation is a well studied problem in mechanics. Currently this problem is addressed with specialized algorithms which can be divided into two classes: refinement approaches and coupled approaches. Refinement approaches alter the number of degrees of freedom (DOFs) in the simulation as a function of performance and accuracy. On the other hand, coupled approaches avoid refinement by using two (or more) simulation techniques which operate on different scales. These methods are common for turbulence simulation in fluids wherein a coarse fluid simulation is coupled (one-way or two-way) to a fine scale turbulence model. However, both refinement approaches and coupling approaches can be difficult for an artist to control as they can affect global changes in the simulation, requiring further parameter tuning to recapture underlying coarse motions. Unlike these previous approaches, what we seek is simulation paintbrush that an artist can use to augment kinematic or coarse dynamic animations with new details in an intuitive fashion.

In this paper we propose a new coupling approach for multi-scale, artist controllable simulation: the Multifarious Hierarchy of Mechanical Models (MHMMs). Multifarious Hierarchies are a general approach that allows a motion to be decomposed using a set of arbitrary, overlapping degrees-of-freedom. Analogous to the layers of popular image editing tools, MHMMs allow an artist to layer spatially varying simulations in order to carefully add secondary detail to coarse motions (Fig. 1). Using a novel kinetic filtering approach, MHMMs avoid kinematic redundancy in the displacement discretization allowing for intuitive augmentation of the motion.

Furthermore, the geometry in the overlapping region experiences a hierarchically combined deformation field and this eases modeling since no geometrical operations are required to achieve the coupling. Finally, MHMMs increase computational efficiency because the total number of system degrees-of-freedom can be optimally chosen by the artist -no more than necessary will be used.

MHMMs allow an artist to produce robust, high-performance, multi-scale simulation results in an intuitive, layered fashionyielding whole simulation results that are greater than a sum of their parts.

Contributions

The major contribution of this work is a general formulation for the layered addition of secondary motion to coarse animations. MH-MMs allow for efficient, detailed simulations and in comparison to other methods, have the following advantages:

• Controllable: Users can add additional detail to a simulation at specified locations and at specified scales. • Spatially-Varying Layers: Deformation models can be applied globally or locally allowing detail to be limited to appropriate regions of the domain (i.e. at contact points etc...). • Arbitrary Number of Levels: Any number of overlapping layers can be used to craft the desired animation. • Layers of Varying Dimension and Material: Combine rods, shells and volumetric models seamlessly; the detail can be used to represent a different material from the coarse model • Efficient: MHMMs can be decomposed into independent systems that can be solved in parallel at runtime • General: Any combination of degrees-of-freedom (particles, rigid and affine frames, modes...) or deformation models can be combined into a hierarchy • Fully-Coupled: Forces and motions applied at one hierarchy level are propagated to the rest of the hierarchy

The remainder of our paper is organized as follows. Background is provided in sections 2 and 3, while our contributions are explained in sections 4 through 5. Section 6 is reserved for implementation details while results are presented in section 7.

Related work

Adding relevant fine-scale details has long been an important area of research in graphics and engineering. As in the introduction we can divide previous methods into refinement approaches and coupling approaches. We begin by reviewing refinement approaches which can themselves be categorized by their degree of hadaptivity and p-adaptivity. H-and p-adaptivity describe the ways in which simulation methods attempt to balance fidelity and performance. H-adaptive methods concern themselves with introducing new degrees-of-freedom by splitting individual elements while p-adaptive methods increase the degree of the polynomial approximation used for field variables. Both approaches have been well studied. In terms of h-adaptivity, [START_REF] Debunne | Interactive multiresolution animation of deformable models[END_REF] introduce a multi-resolution, particle-based methodology for simulating linearly elastic objects. They provide a formalism for dynamically adding and removing particles in order to maintain a desired level of accuracy, and to locally switch between resolutions. Other hadaptive remeshing approaches for Finite Element Simulation have been proposed, most recently in the context of thin shells, elastoplastic simulation and mixed fluid solid simulations [START_REF] Wicke | Dynamic local remeshing for elastoplastic simulation[END_REF][START_REF] Narain | Adaptive anisotropic remeshing for cloth simulation[END_REF][START_REF] Clausen | Simulating liquids and solid-liquid interactions with lagrangian meshes[END_REF][START_REF] Narain | Folding and crumpling adaptive sheets[END_REF]]. However, these techniques are limited to triangular and tetrahedral meshes which are unsuitable for many applications [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]]. [START_REF] Tournier | Velocity-based Adaptivity of Deformable Models[END_REF] introduced a h-adaptive approach for frame-based deformable solids allowing coarsening and refinement. Grinspun et al. [2002] present a h-hierarchical approach which reformulates remeshing as basis refinement. Their method is generally applicable to all element types but is limited to a single mechanical model and aligned levels of detail. This can require complicated mesh coarsening as a preprocessing stage. Padaptive methods have been used with success to simulate fluids on coarse grids [START_REF] Edwards | Detailed water with coarse grids: Combining surface meshes and adaptive discontinuous galerkin[END_REF] however for complicated domains such as solid objects, difficulties can arise. Because the polynomial basis is refined, p-adaptive methods begin with a coarse mesh that reasonably captures the geometry. As such, mesh generation can be problematic [Szabó et al. 2004].

Next we review coupling approaches which attempt to add detail to a simulation by "gluing" differing representations of an object together. These methods can be further divided into surface embedding methods, attachment methods and overlaying methods. Surface methods augment coarse volumetric simulations with highresolution surfaces. Attachment methods provide a mechanism for coupling two different, non-overlapping simulations and overlaying methods use the superimposition of overlapping discretizations to add detail.

Surface embedding methods are popular in graphics and are used to retain reasonable surface detail for rendering while achieving fast simulation speed. Müller et al [2002] use barycentric coordinates to map a coarse deformation to a high resolution surface mesh. Similar approaches have been used to simulate viscoelastic objects with thin features [START_REF] Wojtan | Fast viscoelastic behavior with thin features[END_REF] as well as liquids [START_REF] Wojtan | Liquid simulation with mesh-based surface tracking[END_REF]]. [START_REF] Wojtan | Fast viscoelastic behavior with thin features[END_REF]2011] use the surface to track features but limit its influence on the underlying simulation (limited to computing surface tension forces). Some recent approaches procedurally deform the embedded surface, either from geometric analysis [START_REF] Rohmer | Animation wrinkling: Augmenting coarse cloth simulations with realistic-looking wrinkles[END_REF] or based on examples [START_REF] Wang | Example-based wrinkle synthesis for clothing animation[END_REF][START_REF] Seiler | Enriching coarse interactive elastic objects with high-resolution data-driven deformations[END_REF][START_REF] Zurdo | Animating wrinkles by example on non-skinned cloth[END_REF]. In these methods the fine detail cannot fully respond to user manipulation, relegating it to a more cosmetic role. In general these methods only deal with high-resolution detail at the surface of an object.

Attachment methods have also received a great deal of focus since [START_REF] Sifakis | Hybrid simulation of deformable solids[END_REF] proposed a general methodology to connect models on mesh boundaries. Twigg et al proposed Point Cloud Glue [2010] to easily attach any sets of points using the Procrustes transform. Because simulations of differing resolutions can be attached, this method could be used to resolve details at certain parts of the model. However, attachment methods lack a mechanism for dealing with overlapping discretizations and so dynamic level-ofdetail in a particular simulation domain would require an additional refinement algorithm. Furthermore, there are complications when allowing the detail level to be fully dynamic. Other methods based on substructuring are also commonly applied in engineering and graphics [Barbič and Zhao 2011;[START_REF] Kim | Physics-based character skinning using multidomain subspace deformations[END_REF]. Often per-component reduced models are used to improve performance. Again, an additional algorithm for controlling level-of-detail within each substructure is required.

Overlapping methods have also been explored since the early work of [START_REF] Faloutsos | Dynamic free-form deformations for animation synthesis[END_REF] where local Freeform deformation (FFD) lattices are embedded into larger global ones to allow for fine grained animator control. This method illustrates a particular instance of an MHMM and in our paper we focus on generalizing the concept to layerings of arbitrary discretizations. Two-way coupled algorithms were initially explored by Terzopoulos and collaborators [1988;[START_REF] Terzopoulos | Dynamic 3d models with local and global deformations: deformable superquadrics[END_REF] ing position constraints specific to the behavior of the stiff hull. [START_REF] Harmon | Subspace integration with local deformations[END_REF] enrich a space of modal deformations (see Barbic et al [2005]) with analytically defined detail functions (for a point force applied to an elastic half-space) to increase the expressivity of the model. These methods are very much in line with the spirit of our work but are, again, limited to particular instances of MHMMs (rigid modes and deformable superquadrics, poking functions adding details over modal models). Finally, the Eulerian-on-Lagrangian method [START_REF] Fan | Eulerian-on-lagrangian simulation[END_REF][START_REF] Fan | Active volumetric musculoskeletal systems[END_REF] allows the coupling of a fine scale Eulerian simulation to any coarse scale Lagrangian simulation. The method is general but limited to Eulerian-Lagrangian coupling. Additionally the Eulerian simulation must cover entire deformable objects, meaning that detail cannot be added locally thus negating any potential performance gains (i.e in [START_REF] Fan | Active volumetric musculoskeletal systems[END_REF] entire groups of muscles are covered with simulation grids rather than one muscle being covered by many grids) . Conversely, the focus of MHMMs is to provide a principled, controllable, layered methodology for adding local detail to coarser simulations, something that none of the previous work in this area provide.

Embeddings

This section provides background on embeddings, introduces notations and motivates our approach. Embeddings are intuitive -if one is playing a game of darts, the darts themselves become embedded in the board after each toss. Moving the board then causes all the darts to move with the same rigid motion.

A mathematical embedding (often used for animation [Capell et al. 2002;[START_REF] Galoppo | Fast simulation of deformable models in contact using dynamic deformation textures[END_REF]) formalizes this idea that the motion of an embedded mathematical object (a 3D point, deformation gradient, etc..) is influenced by the movement of some encapsulating domain (such as a rigid transform, a per-vertex displacement field or an affine blending). More generally, any mapping from one set of DOFs to another, such as vibration modes to local strains, can be used to build an embedding. Embeddings are particularly useful for synchronizing the motion of the embedded object with that of the encapsulating one. In this paper we use the term master to refer to the DOFs of the encapsulating object (the dart board) and the term slave to refer to the DOFs of the embedded object (the darts). Slave DOFs are kinematically "slaved" to master DOFs, hence their name. Furthermore we use the term coarse to refer to master DOFs that parameterize a lower-resolution simulation layer while the terms fine or detail refer to master DOFs that parameterize higher-resolution, detail layers.

Let us consider a set of material points, with X the vector of all their material coordinates, constant over time. As the state vector of the object, q (containing the displacements of the master DOFs from their initial coordinates), changes over time t, so do the displacements of the embedded, or slave, points, ue. The mapping relation is:

ue(t) = J (q(t), X) (1)
For the sake of clarity, we will drop the dependence on t in subsequent notations. Slave DOFs may in turn be the masters of other sets of DOFs, forming a hierarchy with the root as only set of master DOFs (Fig. 2). The velocities ve at the embedded level are given by: ve = J q (2)

where J = @J @q is the Jacobian matrix of the mapping. Note that positions and velocities are only propagated "top-down", from the master to the slave DOFs. The Jacobian matrix is generally rectangular, and so mapping the velocities in the opposite direction requires solving an inverse kinematics problem, which often has no unique solution, if any. Conversely, the forces and impulses are only propagated bottom-up. Let fe be forces at the embedded level, and fq the corresponding forces at the master level. The Principle of Virtual Work states that the power of the forces must be the same at the two levels: v T e fe = qT fq, for any value of ve. This, combined with Eq.2, implies that fq = J T fe.

(3)

Mass and stiffness matrices can be transferred "bottom-up" (from slave to master DOFs). For instance, if the inertia of the object is defined as a matrix Me, the equivalent matrix at the top level is Mq = J T MeJ. While the complete derivation of all the properties of embeddings is out of the scope of this paper, we note that the above relations between forces and displacements at different levels of the kinematic hierarchy hold for any set of generalized coordinates and forces.

Hierarchies of Embeddings

Principle

MHMMs consist of augmented embedding hierarchies. Rather than maintaining a strict master-slave relationship between embedded DOFs, we instead augment the embedding with its own master DOFs so that it can contribute new motion to the animation. One way to perform this augmentation would be to use a simple sum of embeddings ue(t) = J1(q1(t), X) + J2(q2(t), X). For small displacements this produces acceptable results, but its uncorrelated terms result in unacceptable artifacts as soon as large deformations or rotations occur.

In order to build an augmented embedding suitable for animating large deformations we instead modify Eq.1 by using a second embedding to change the material coordinates used by the first:

ue = Jc(qc, X + J d (q d , X)) (4)
where Jc(qc, X) is a coarse displacement field and J d (q d , X) is a new model used to add detail, as if a deformable object were embedded in the coarse model. In this approach, the detail changes the material coordinates used in the coarse model. We stress that such a hierarchy of embeddings can be superimposed on any mesh in order to parameterize its motion in a layered fashion.

Fig. 3 shows the two-level case of an MHMM while Fig. 4a illustrates the relationship between embeddings, DOFs and relevant physical objects (such as quadrature points). The top row of Fig. 4a shows the relationship between the coarse DOFs (c) and a detail layer consisting of a local FEM simulation (d). All other objects One potential problem with this approach is that combining embeddings using Eq.4 may require intense computations at each time step to update the coarse embedding based on the modified material coordinates, and may even introduce discontinuities when an embedded point moves from one cell of a mesh-based coarse model to another. To alleviate this problem, we express the detail motion in a local frame associated with each embedded point. This is essentially a first-order expansion of Eq.4:

ue = Jc(qc, X) + FcJ d (q d , X) (5)
where Fc is a block-diagonal matrix. For each embedded point i, the diagonal block represents the deformation gradient

F i c = @J i c
@X i + I of the coarse displacement at embedded point i. This non-orthonormal basis represents the rotation and deformation created by the coarse displacement, in which the detail displacement is applied, as illustrated by the example in Fig. 3. Here, we assume that the coarse deformation field is locally affine (i.e.,Fc is uniform around X), which is exact for some coarse models such as triangular or tetrahedral finite elements. In the general case, this assumption is reasonable when detail displacements are small or when coarse deformations are smooth. In contrast, summing embeddings requires small displacements of all models. Following Eq.2, the embedded velocity can be expressed as:

ve = J q = ⇥ Jc J d ⇤  qc qd (6)
where the Jacobian is a sparse matrix with only two non-null blocks per embedded point:

J i c = [ @J i c @qc + @F i c @qc J i d ],
and

J i d = [F i c @J i d @q d ]
for the coarse and detail models, respectively. However, it is easy to see that straightforwardly combining arbitrary embeddings using this method is not physically correct. Naïvely propagating a force from the embedded layer to both the coarse and detail layers (Eq.3) results in its duplication. Forces must be properly partitioned and then applied to each layer, ensuring that the total sum of work performed on the object is consistent.

Similar issues arise when dealing with other intrinsic object properties. Properties such as mass and stiffness are defined independently of the kinematic discretization either analytically or on a fine voxel map. These properties are sampled using embedded particles, as illustrated in Fig. 4a. A given mass matrix Me at the embedded level is not directly exploitable in Newton's law f = Ma, since the latter is only valid for independent DOFs, which requires us to modify the mass matrix at a particular hiearchy level. Using the Jacobian of Eq.6, the equivalent mass matrix at parent level would be:

M = J T MeJ =  J T c MeJc J T c MeJ d J T d MeJc J T d MeJ d (7)
Fan et al. [2013] observed that, while the diagonal blocks are symmetric-positive-definite (SPD), this mass matrix is generally singular, due to the off-diagonal coupling blocks. This happens in the case of redundancies between the coarse and detail embeddings.

For instance, if the coarse field is a rigid displacement and the detail is an FEM mesh, then a translation of the embedded object can be obtained using a translation of the frame as well as a uniform displacement of the FEM nodes. To alleviate the singularity, they use a least-squares approach to transfer fine scale Eulerian displacements onto a coarse Lagrangian discretization. In contrast to this, we propose a more physically-based approach, using the mass matrix, as explained in the next section.

Kinetic Filtering

We attempt to find the corresponding generalized velocities at the coarse level, qc, that best match ve, the velocity at the embedded level. This can be done by solving a mass-weighted least squares problem, the solution of which yields qc = M 1 c J T c Meve, where Mc = J T c MeJc. Note that, if one takes Me to be a multiple of identity, this reduces to the classical pseudo inverse solution to a system of equations. Now we can define v 0 e the embedded velocity captured by qc as

v c e = Pve (8) P = JcM 1 c J T c Me. (9)
We can also compute the remaining detail level velocity as where I is the identity matrix. This leads us to define the complementary operator P = I P which computes detail velocity (at the embedded level) from the total velocity. It is straightforward to see that p c e = Me Pve yields an embedded momentum in the nullspace of J T c . Concretely, if the coarse layer is a rigid transform and the fine layer is a finite element simulation, P strips rigid motion out of the FEM velocity field while maintaining the non-rigid motions. This avoids kinematic redundancy since rigid motions can only be manifested at the coarse layer.

In order to continue, we define the Jacobian of our system as

J = ⇥ Jc PJ d ⇤ (11)
and consequently the momentum computed by the off-diagonal blocks of the generalized mass matrix for an arbitrary detail velocity q⇤ d is

J T c Me PJ d q⇤ d = J T c p ⇤ e = 0. (12)
Because q⇤ d was chosen arbitrarily, this is equivalent to the offdiagonal mass matrix block, J T c Me PJ d , being 0. Note that this does not occur when using the non-mass-weighted formulation of [START_REF] Fan | Eulerian-on-lagrangian simulation[END_REF] and that later, we will exploit this property to evaluate hierarchy levels in parallel.

The resulting form of Newton's equation is:

 J T c MeJc J T d PT Me PJ d  ac a d =  J T c fe J T d PT fe (13)
The upper-left block is SPD provided that Me is symmetric and Jc is not rank-deficient. While the lower-right block is singular due to the filter, it can be solved using a Filtered Conjugate Method [START_REF] Ascher | On the modified conjugate gradient method in cloth simulation[END_REF]].

An example of tetrahedral FEM detail combined with a rigid coarse motion is shown in Fig. 5. The rigid frame is fixed for illustration purposes, so that only detail deformations are shown. An external force applied to a single point creates deformations. The force is split by the transposed Jacobian, and the FEM model gets a force distribution which can only generate null linear and angular momenta. It can thus generate any displacement but a rigid one. Notice that there is no fixed point in our detail model. Contrary to previous approaches, the connection between the two models is made by the filtered Jacobian, not by attachment points. This allows us to straightforwardly combine linear modal deformation modes as detail on top of rigid motion, as shown in Section 7.4, to avoid the computation of nonlinear deformation modes.

Forces

Our framework is compatible with all the usual material laws and contact behaviors. In the simple approach illustrated in Fig. 4a, the surface layer s, as well as the deformation energy integration samples (Gauss points in layer v), are embedded in layer e using interpolation. Embedding deformation gradients in layer v allows us to implement all the standard strain measures (such as corotational and Green-Lagrange), and associated constitutive laws (such as Hookean or hyperelastic). The associated generalized forces are stresses, which are propagated as forces to layer e, which also collects the contact forces mapped from p through s. These forces are split between coarse and detail forces as explained in Section 4.2.

N-levels hierarchies

Kinetic filtering can be straightforwardly extended to handle more than two levels in the hierarchy. Let qi, Ji(t)(qi, X) and Fi 1 = @J i 1 @X i 1 respectively the coordinates of the DOFs, the detail and the local frame at level i. Each detail layer is applied in the parent frame, so that for N levels, we get:

ue = J1 + F1J2 + F1F2J3 + . . . = J1 + P N k=2 ( Q k 1 j=1 Fj)J k (14)
The corresponding Jacobian is:

J = ⇥ J1 . . . Pi 1Ji . . . ⇤ (15) P0 = I, Pi = Pi 1.(I Pi) (16) Pi = Pi 1JiM 1 ii J T i PT i 1 Me Pi 1 (17) Mii = J T i PT i 1 Me Pi 1Ji (18)
This formulation achieves decoupling since off-diagonal blocks of the mass matrix are null (see proof in Appendix A).

Optimizing Performance

Constructing the kinetic filter requires relating momenta of one level of the hierarchy to the momenta of another. In the naïve implementation described above we do this by storing high resolution information on the embedded particles in layer e (the density of which must be sufficient to represent the deformation field at the most detailed layer, d, in the range of the detail deformation). Thus the matrix P, used by the kinetic filter, is dense and is of dimension of the number of particles in the embedded layer. While P is never inverted, this may result in increased computation times due to the cost of performing dense matrix-matrix multiplication to form the per-layer mass matrices (the size of these is strictly determined by the number of master DOFs in the particular layer). In order to reduce the computation time without modifying the sampling scheme, we present an optimization to bypass the embedded layer e completely, thus directly connecting v and s to the master layer c-d. Each of the points in layers s and v have their own material coordinates and can thus be directly embedded into c-d.

For layer v, we compute deformation gradients Fv by spatially differentiating the embedding of Eq.5 with respect to the material coordinates of the sampling points:

Fv = FcF d (19)
This corresponds to the composition of the coarse and detail deformation gradients as shown in Fig. 3. The stresses fv are computed using the deformation gradients and the local material constitutive law, and integrated in space using standard cubature. Stresses could be directly mapped up as forces on layers c and d based on the transposed Jacobian matrix (Eq.3). However, the question is how to filter the Jacobian of this embedding, to avoid "double counting" stresses while converting them to forces at the master level.

Let Jcv be the Jacobian of the embedding of deformation gradients in the coarse model, and J dv be the unfiltered Jacobian in the detail model, obtained by differentiation of Eq.19 with respect to the coordinates of the DOFs. We convert stresses to filtered forces in the following way:

f d = J T d PT Me[JcM 1 c J T cv + PJ d M 1 d J T dv ]fv (20)
The above equation transforms the force fv from v to d using the following steps (right-to-left):

1. Transform forces from v to c and d using J T cv and J T dv 2. Convert forces at d and c to accelerations using M 1 d and M 1 c 3. Sum contributions to the corresponding embedded accelerations using Jc and J d 4. Compute the force at the embedded level using Me 5. Split the force on the master layer using kinetic filtering By re-arranging the above equation using Eq.9, we get:

f d = (I P 0T )J T dv fv (21)
P 0 = M 1 d J T d MePJ d (22) J 0 = ⇥ Jcv J dv (I P 0 ) ⇤ ( 23
)
Note that now the dimension of our kinetic filter, P 0 are dictated by the number of DOFs in the detail level only, a signfigant reduction in size that yields a corresponding performance increase

The square matrix (I P 0T ) acts as a filter on the detail forces, while its transpose acts as a filter on the detail velocities. Moreover, since the above formulation holds for any other slave layers such as s, it can be applied once to the force on d and then accumulated bottom up, as illustrated using d 0 in Fig. 4b.

Decoupled implicit integration

For implicit time integrators we can gain further performance improvements by modifying our kinetic filter so that it decouples hierarchy levels during the implicit solve. This allows us to parallelize our velocity update steps on a per-layer basis. Consider the Implicit Euler method, which requires repeated solution of the following equation:

M h 2 K a = f + hKv ( 24
)
where K = @f @q is the stiffness matrix, and h is the time step. Given the mass matrix Me at the embedded level, the mass matrix at the master level is computed as explained in Section 4.2 using the filtered Jacobian of Eq.11. We do not assemble the stiffness matrix at the master level, since we use a Conjugate Gradient solver and propagations to compute the matrix-vector products. Let H = M h 2 K be the integration matrix. The filtering method to cancel the off-diagonal blocks of the mass matrix presented in Section 4.2 can straightforwardly be extended to this matrix by replacing the projection P of Eq.8 with:

v 00 = Sv = JcH 1 c J T c Hev, (25)
Equation 24 becomes:

 J T c HeJc J T d (I S T )He(I S)J d  ac a d =  J T c (fe + hKev) J T d (I S T )(fe + hKev) (26)
This results in two systems coupled only through the right-hand term. These can be solved independently, possibly using specialized solvers or in parallel. Our coupling is neither fully explicit nor fully implicit. The matrix S is used to implicitly account for the change in coupling force due to changes in the detail displacements. However, it ignores changes in the kinetic filter and J d . This is analogous to standard implementations of co-rotational FEM, which assume that the rotation matrix is fixed across time steps. In practice we've found this integration scheme to be stable for large time steps. An example is shown in Section 7.4 wherein the rigid dynamics are handled by a generic solver while the modal model is handled using a specialized, preconditioned conjugate gradient solver described in the following sections.

Pre-factored detail systems

In MHMMs, detail motions are expressed in local frames relative to the coarse motion (J T d in the right hand side of eq 26 transforms filtered forces to the detail coordinate system, where the detail system is solved). The nature of the detail motion is that it is highfrequency, typically with lower amplitude -thus we can apply small deformation simplifications. By pre-factoring the lower right term of the matrix of Eq. 26, supposing that J d is small enough, we can increase computational time, at the cost of accuracy. Examples of such systems are shown using linear FEM(Figure 7), and linear modal simulations(Figure 10).

For non-interactive scenarios (e.g. animations for films) this prefactoring has no drawbacks. However, in applications wherein the user manipulates the mesh by applying an external force, the prefactorization prevents us from implicitly integrating this interaction force. It is worth nothing that this limitation is common to most methods that prefactor the implicit integration system matrix.

Preconditioning Nonlinear Problems

Our prefactored, per-layer, system matrices can also be used to accelerate the solution of non-linear detail simulations by preconditionning the detail solver using the initial prefactored system matrix (see Figure 8): Performance improvements of up to 50% have been observed in our experiments. This greatly reduces the cost of adding extra detail to any simulation.

Precomputation Cost

Our precomputation requires the assembly of projection matrices (involving the inversion of mass matrices as in standard direct solvers, and a matrix multiplication of the size of the embedding at most) and the proposed prefactorization (cholesky factorization on a matrix of the layer size). Practically speaking, precomputation time is on the order of a few simulation timesteps at most.

Implementation Details

In this section we review some relevant implementation details as well as expound upon important optimizations that can be applied to MHMMs.

At initialization time, we first organize the model based on the property map of the object, as illustrated in Fig. 4a. Computing the mass Me at the embedded level and storing it as a constant in the master layer is a reasonable and efficient approximation. Based on the stiffness distribution and the embedding, we compute the volume samples of layer v. We then set up the direct embeddings and remove the embedded layer e as explained in Section 5 and illustrated in Fig. 4b. Since gravity generates translational forces, its contribution to the detail is null when translation is kinematically feasible by the coarse model. In this case, we apply it directly at the coarse level. Otherwise, it is possible to split the weight across volume sampling points.

In general, the Jacobians change over time and must be updated at each step of the simulation, and the filters as well. However, many popular blendings such a linear FEM, linear modes, or affine frame-based blendings have constant Jacobians, and result in constant @Jc/@qc or @J d /@q d matrices which allow us to save computation time.

Assuming that the fine displacement of the detail J d is small, the coarse Jacobian matrix Jc and filter P can be considered constant. When the coarse displacement is large, the change of local frames Fc in Eq.5 requires updates of J d .

Lazy updates of Jacobians can lead to poor filtering performance and this manifests as coarse motion bleeding into the detail layers. This is due to violating the constraint that off-diagonal blocks of the mass matrix should be zero. In practice we correct this using standard stabilization techniques. Such stabilization may be achieved by fixing the edges of the detail model, which also ensures continuity in the displacement, at the interface of this model. With stabilization in place the method is robust enough so that, in practice, any linear solver like CG or MINRES will work for every scenario.

Pre-computing the filter once for the whole simulation does not lead to visual artifacts (as can be seen in the accompanying video) while recomputing the filter doubles the computation time of a given time step.

Results

In this section, we exhibit the versatility of our approach using different combinations of deformation models (2d/3d finite elements, rigid, frame-based, modal subspace) in conjunction with artist positioned secondary motion.

Physical Plausibility

To show that our algorithm produces physically plausible results we compared it to a high-resolution FEM simulation of a beam undergoing 3-point bending. Figure 6 shows both the FEM result and the result of an MHMM composed of two coarse hexahedra, manually overlaid with a fine hexahedral grid. The relative displacement error is less than 20% and visually the results are indistinguishable.

FEM on a rigid frame

An artist can choose to apply detail motion to a whole object, as shown in Fig. 7. Here, a rigid model is enriched with tetrahedral FEM (red edges on Fig. 7-left), covering the entire object. Note that there is no attachment point between the rigid and the FEM models. Even for a very stiff material, the weight is correctly simulated in these two examples, because the kinetic filter maps it to the rigid model. This alleviates the well-known gravity artifact which occurs when iterative solvers are used with stiff FEM. Incomplete solutions lead to stiff objects falling more slowly.

FEM on frame-based deformation models

Our method allows an artist to add local secondary motion to any type of coarse simulation or kinematic animation. In Fig. 1 the coarse model of the octopus is frame-based [START_REF] Gilles | Frame-based Elastic Models[END_REF], which allows for efficient simulation of the global motion using linear blend skinning. The artist then chooses to add a small, highresolution FEM patch to the octopus' head in order to capture the deformation caused by its interaction with the hook.

Figure 9 shows a similar example where the artist has added a cylinder of high-resolution FEM to capture the pinching caused by the contacting string. Allowing artists to easily had detail where its needed has the added benefit of reducing computational time as well as offering a corresponding reduction in memory footprint (Table 1).

Non-linear Modal Subspace

Non-linear subspace deformation can be easily and efficiently achieved with our technique by combining a coarse frame-based deformation field and linear subspace deformation modes. These modes can be drawn from simulation data, capture data or hand crafted by an artist. The MHMM locally transforms each mode by the blended displacement of the frames. The well known artifacts due to linearity are reduced as the number of frames increases. This is akin to a co-rotational simulation method wherein the rotation is encoded by the frames rather than resulting from a polar decomposition.

Again we see that MHMMs allow for better control of animations. As shown in the accompanying video, the scripted displacement of one of the frame can be automatically enriched with simulated deformations (based on the remaining frames and modal degrees of freedom).

To exhibit particular visual features, a third level can be added. In Figure 10, local detailed deformations are achieved using a local FEM patch. In Fig. 11, we augment a global, linear modal model using local subspace deformation modes which allows for fine tuning of the deformation degrees of freedom.

Multi-dimensional material

MHMMs also provide a way for artists to modify the material composition of simulated object. For instance, because our approach makes no assumption on the topology of the models, one can easily add a stiff skin to coarse soft tissue (Figure 12)

The coarse FE achieves volume preservation, while the fine surface creates folds. The coupling achieved by the kinetic filter generates the well-known buckling behavior without additional constraints. In Fig. 13, we simulate the pinching of a patch of skin with hypoderm using the epiderm as 2D detail on top of a two-material coarse volumetric grid.

Multilayer Editing of Animations

Here we show an example of multilayered animation editing, akin to the layered approach used in many popular image editing suites. Our kinetic filter helps ensure that adding new detail levels to a simulation will have a limited effect on any underlying gross motion. Figure 14 shows a 3 level hierarchy applied to an animated character. The coarse motion is prescribed by an artist using standard skinning techniques. The artist has then added a coarse FEM grid to capture the jiggling of the character's belly along with a finer, local grid to add the indent caused by the character's finger. Our approach can emulate Rig-Space physics [START_REF] Hahn | Rig-space physics[END_REF]] using a hierarchy where the master layer is the rig that controls a fine tetrahedral FEM model. However, we improve upon this technique by allowing artists to add physical details without the need to modify the characfer rig itself. One could use a full high resolution simulation such as in [START_REF] Mcadams | Efficient elasticity for character skinning with contact and collisions[END_REF], however using MH-MMs requires far fewer degrees-of-freedom and allows for direct augmentation of the skinned surface animation (rather than being limited to having the bone motion drive a soft elastic character). All while permitting significantly more user control of deformation placement. Not to mention that our technique is compatible with a wider range of material models.

Performance

Frame rates for our examples are given in Table 1. For a fair comparison, we use the same solvers and contact handling mechanisms for all deformable models. MHMMs achieve large speedups compared to models with uniformly fine resolutions (corresponding to the detail we wish to simulate). As is the case with other techniques, the number of degrees of freedom impacts the performance of our method but, we gain invariance with respect to the spatial dimension of the detail. Therefore, for a given number of degrees of freedom, details at arbitrarily fine scales can be simulated in constant time. Other methods would require extremely dense discretizations, re-meshing or adaptive refinement to accomplish this.

Limitations and Future Work

The main bottleneck of MHMMs is potential coupling of all detail DOFs through the coarse layer. This results in a dense kinetic filter matrix, which slows down the required matrix products. In such cases, a more standard, high resolution simulator is likely to achieve better performance. For MHMMs, pre-factorization mitigates this problem, as the filter is only required to compute the right hand side vector. We suspect that filtering may be further accelerated using a multigrid approach, but defer this study to future work. Finally, the mechanical coupling between layers could be improved using the geometric stiffness [START_REF] Tournier | Stable Constrained Dynamics[END_REF]] that accounts for the variation of the Jacobian within timesteps.

Conclusion

Multifarious Hierarchical Mechanical Models are a new approach to solid simulation, based on hierarchically coupling deformable models. Using a novel operator, the kinetic filter we provide a single tool that can be used by artists to allow simple editing of a simulation as well as to simulate heterogeneous materials (see our skin example) or provide additional detail in standard simulations, even ones which are highly non-linear. As such, MHMMs bring the convenience of layered editing, common place in the image processing domain, to physically-based animation.

Our algorithm is implemented using the SOFA open-source library supported by INRIA and partners [START_REF] Faure | SOFA: A Multi-Model Framework for Interactive Physical Simulation[END_REF], and is publicly available in the MultifariousHierarchies plugin (https: //gforge.inria.fr/projects/mhmm).

Scene Octopus (Fig. 1) Walrus (Fig. 10) Tourniquet (Fig. 9)

DOFs

A Off-diagonal blocks for N-levels

In the expression of an upper off-diagonal block Mij, i < j, we split Pj 1:

Mij = J T i PT i 1 Me Pj 1Jj = J T i PT i 1 Me Pi 1(I Pi) j 1 Y k=i+1 (I P k ) ! Jj
The first part of this formula is null, due to the properties of our filters:

J T i PT i 1 Me Pi 1 (I Pi) = J T i PT i 1 Me Pi 1 J T i PT i 1 Me Pi 1 Pi 1 Ji | {z } M i M 1 i J T i PT i 1 Me Pi 1 = 0
The property Pi 1 Pi 1 = Pi 1 holds since Pi 1 is a projection.

  (a) Rest state (b) FEM deformation onlyF c (c) Frame-based deformation only

Figure 3 :

 3 Figure 3: Example of a two-level MHMM, using a frame-based coarse model (blue) and an FEM detail model (red). The detail displacements J d (3b) are applied in local frames defined by the deformation gradients Fc of the coarse displacements (3c) to generate the combined deformations (3d).

Figure 4 :

 4 Figure 4: Structure of the two-level MHMM shown in Fig.3.

Figure 5 :

 5 Figure 5: Filtered detail motion. Left: rest shape. Middle, right: detail deformations generated by external forces.

Figure 6 :

 6 Figure 6: Comparison on flexion. Left: standard FEM. Right: MHMM with coarse FEM + detail In the accompanying video, the simulation of a cantilever beam in contact shows a major advantage of MHMMs with regards to artist directed level-of-detail. By adding detail only in desired locations (near the contact point) MHMMs keep the global motion of the beam close to the input coarse motion. The fine FEM grid leads to spurious secondary motion that alters the beams trajectory signifigantly. Note that all the beams reach identical final configurations.

Figure 7 :

 7 Figure 7: Rigid frame (blue) + global FEM (red). Left : initial position. Middle : simulated models shown separately. Right : Combined result. Grey objects are fixed colliders.

Figure 8

 8 Figure8shows an example of a local FEM volumetric patch hierarchically attached to a moving frame. This simple simulation, which could also be produced using traditional coupled models, shows that our detail model can undergo large deformations. The detail model is fixed at its boundary with the rest of the object, to enforce the continuity of the displacements.

Figure 8 :

 8 Figure 8: Rigid frame + local FEM: a baby bottle with a deformable soother. The detail level FEM is displayed with its local coordinate system centered on the origin.

Figure 9 :

 9 Figure 9: A string strangulates a deformable cylinder (purple). A local tetrahedral FEM patch (red) allows a precise deformation in the contact area while coarse affine frames (blue) simulate the global deformation.

Figure 10 :

 10 Figure 10: A 3-levels hierarchy combining frame-based, linear modal subspace, and hexahedral finite elements models. Left: Frame and FEM discretization; Middle: volume samples; Right: Simulation.

Figure 11 :

 11 Figure 11: Frame-based, global subspace and local subspace deformation models combined in a 3-level hierarchy.

Figure 12 :

 12 Figure 12: A detailed thin shell embedded in a single coarse hexahedral element.

Figure 13 :

 13 Figure13: Skin folding simulation using a coarse volumetric FEM mesh (in blue) combined with a thin shell patch and several material properties.

Figure 14 :

 14 Figure 14: A full body with a coarse rigid motion and three levels of hierarchy easily controlled by an artist.

Table 1 :

 1 10 affine 6878 points 10 affine + 48 points 4566 points 2 rigids + 16 modes 2 rigids + 16 modes + 111 points 7500 points 5 affine + 510 points Performances GRINSPUN, E., KRYSL, P., AND SCHR ÖDER, P. 2002. Charms: a simple framework for adaptive simulation. In ACM Transactions on Graphics (Proc. SIGGRAPH), vol. 21.

	Integration Points	10	37912	154	26392	50	390	36450	2440
	Frame Rate (Hz)	60	4	35 (pre-factorized)	1	100	15 (pre-factorized)	<0.1	2 (pre-factorized)
				25 (full solve)			10 (full solve)		2 (full solve)

Acknowledgments

The authors wish to thank Laura Paiardini for her help with the models, textures, renderings and video editing and Romain Testylier for his help with the rod simulation. We also thank the anonymous reviewers for the suggestions for improving the paper. This work was supported by funding from French FUIs Dynamit and Collodi, French ANR project SoHusim and INRIA project Morphogenetics.