Jean-Michel Muller
email: jean-michel.muller@ens-lyon.fr

Scaling Newton-Raphson division iterations to avoid double rounding

When performing divisions using Newton-Raphson (or similar) iterations on a processor with a floating-point fused multiply-add instruction, one must sometimes scale the iterations, to avoid over/underflow and/or loss of accuracy. This may lead to double-roundings, resulting in output values that may not be correctly rounded when the quotient is in the subnormal range. We show how to avoid this problem.

introduction

Throughout the paper, we assume a radix-2, precision-p, floating-point system that is compliant with the IEEE 754-2008 Standard for Floating-Point Arithmetic [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF]. We also assume that a fused multiply-add (FMA) instruction is available. That instruction evaluates expressions of the form xy + z with one final rounding only. We also assume that the ambient rounding mode is round to nearest (this is the only one for which the problem we are dealing with, namely double rounding, occurs).

Many algorithms have been suggested for performing divisions. Here, assuming we wish to evaluate the quotient b/a of two floating-point numbers, we focus on algorithms that first provide an initial approximation q to the quotient and an approximation y to 1/a, and refine it using a "correcting step" [START_REF] Markstein | Computation of elementary functions on the IBM RISC System/6000 processor[END_REF]. r = RN (b -aq), q ′ = RN (q + ry),

where RN(u) is u rounded to nearest (even). Under some conditions made explicit in Theorem 2-roughly speaking, if q and y are close enough to b/a and 1/a, respectively, and no underflow occurs-, then q ′ = RN(b/a).

1

In most applications of that property presented in the literature, the approximations y and q are obtained through variants of the Newton-Raphson iteration, but they might as well result from other means. What matters in this paper is that the correcting step (1) is used.

That correcting step method is applicable only under some conditions. More precisely, what makes the method working is the following lemma, which shows that, still under some conditions, r = b -aq exactly. That lemma can be traced back at least to Markstein's work [START_REF] Markstein | Computation of elementary functions on the IBM RISC System/6000 processor[END_REF]. The presentation we give here is close to that of Boldo and Daumas [START_REF] Boldo | Representable correcting terms for possibly underflowing floating point operations[END_REF][START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF].

Lemma 1 (Computation of division residuals using an FMA) Assume a and b are precision-p, radix-2, floating-point numbers, with a = 0 and |b/a| below the overflow threshold. If q is defined as

• b/a if it is exactly representable;
• one of the two floating-point numbers that surround b/a otherwise; then b -aq is exactly computed using one FMA instruction, with any rounding mode, provided that e a + e q ≥ e min + p -1, and

q = α or |b/a| ≥ α 2 , (2)
where e a and e q are the exponents of a and q and α = 2 emin-p+1 is the smallest positive subnormal number.

For this result to be applicable, we need e a + e q ≥ e min + p -1. This condition will be satisfied if e b ≥ e min + p. Other conditions will be needed for the correcting iterations to work (see Theorem 2 below). Also, the intermediate iterations used for computing q and y may require the absence of over/underflow. All this gives somewhat complex conditions on a and b, that can very roughly be summarized as "the quotient and the residual r must be far enough from the underflow and overflow thresholds". More precisely, Theorem 2 (Peter Markstein [START_REF] Markstein | Computation of elementary functions on the IBM RISC System/6000 processor[END_REF][START_REF] Cornea | Correctness proofs outline for Newton-Raphson-based floating-point divide and square root algorithms[END_REF][START_REF] Markstein | IA-64 and Elementary Functions: Speed and Precision[END_REF][START_REF] Harrison | Formal verification of IA-64 division algorithms[END_REF]) Assume a precision-p binary floatingpoint arithmetic, and let a and b be normal numbers. If

• q is a faithful approximation to b/a, and • q is not in the subnormal range, and

• e b ≥ e min + p, and • y approximates 1/a with a relative error less than 2 -p , and

• the calculations r = •(b -aq), q ′ = •(q + ry)
are performed using a given rounding mode •, taken among round to nearest even, round toward zero, round toward -∞, round toward +∞, then q ′ is exactly •(b/a) (that is, b/a rounded according to the same rounding mode •).

A natural way to make sure that the conditions of Theorem 2 be satisfied is to scale the iterations. This can be done as follows: a quick preliminary checking on the exponents of a and b determines if the conditions of Theorem 2 may not be satisfied, or if there is some risk of over/underflow in the iterations that compute y and q. If this is the case, operand a, or operand b is multiplied by some adequately chosen power of 2, to get new, scaled, operands a * and b * such that the division b * /a * is performed without any problem. An alternate, possibly simpler, solution is to always scale: for instance, we chose a * and b * equal to the significands of a and b, i.e., we momentarily set their exponents to zero. In any case, we assume that we now perform a division b * /a * such that:

• for that "scaled division", the conditions of Theorem 2 are satisfied;

• the exact quotient b/a is equal to 2 σ b * /a * , where σ is an integer straightforwardly deduced from the scaling.

Assuming now that the scaled iterations return a scaled approximate quotient q * and a scaled approximate reciprocal y * , we perform a scaled correcting step r = RN (b * -a * q * ), q ′ = RN (q * + ry * ),

Notice that q ′ is in the normal range (i.e., its absolute value is larger than or equal to 2 emin ): the scaling was partly done in order to make this sure. If 2 σ q ′ is a floating-point number (e.g., if |2 σ q ′ | ≥ 2 emin ), then we clearly should return 2 σ q ′ . The trouble is when 2 σ q ′ falls in the subnormal range: we cannot just return RN(2 σ q ′ ) because a double rounding phenomenon might occur and lead to the delivery of a wrong result. Consider the following example. Assume the floating-point format being considered is binary32 (that format was called single precision in the previous version of IEEE 754: precision p = 24, extremal exponents e min = -126 and e max = 127). Consider the two floating-point input values:

b = 1.00000000001100011001101 2 × 2 -113 = 8394957 10 × 2 -136 , a = 1.00000000000011011001100 2 × 2 23
= 8390348 10 .

The number b/a is equal to 0.1000000000010010000000000000101101001111011001100100000010 This phenomenon may appear each time the scaled result q ′ , once multiplied by 2 σ , is exactly equal to a (subnormal) midpoint, i.e., a value exactly halfway between two consecutive floating-point numbers. Notice that if we just use this scaled result without any other information, it is impossible to deduce if the exact, infinitely precise, result is above or below the midpoint, so it is hopeless to try to return a correctly rounded value. Fortunately, values computed during the last correction iteration will contain enough information to allow for a correctly rounded final result, as we are now going to see.

Avoiding double rounding

As stated in the previous section, we assume we have performed the correcting step: r = RN (b * -a * q * ), q ′ = RN (q * + ry * ), and that the scaled operands a * , b * , as well as the approximate scaled quotient q * and scaled reciprocal y * satisfy the conditions of Theorem 2. We assume that the scaling was such that the exact quotient b/a is equal to 2 σ b * /a * . We assume that we are interested in quotients rounded to the nearest ("even", yet with round to nearest "away" the reasoning is not so different). To simplify the presentation, we assume that a and b (and, therefore, a * , b * , y * , q * and q ′ ) are positive (separately handling the signs of the input operands is straightforward). Since q * is a faithful approximation to b * /a * , we deduce that

q -< b * a * < q + ,
where q -and q + are the floating-point predecessor and successor of q * . Also, since q ′ = RN(b * /a * ), we immediately deduce that q ′ ∈ {q -, q, q + }. This is illustrated by Figure 1.

excluded

q - q * q +
area where b * /a * can be located, the endpoints q -and q + being Figure 1: The number q * is a faithful rounding of b * /a * : this means that q -< b * /a * < q + , where q -and q + are the floating-point predecessor and successor of q * . As stated before, the "double rounding" problem occurs when 2 σ q ′ is a (subnormal) midpoint of the floating-point format. In such a case, to return a correctly rounded quotient, one must know if the exact quotient b/a is strictly below, equal to, or strictly above that midpoint. Of, course, this is equivalent to knowing if b * /a * is strictly below, equal to, or strictly above q ′ . Lemma 1 says that r = b * -a * q * exactly. Therefore, when 2 σ q ′ is a midpoint:

1. if r = 0 then q ′ = q * = b * /a * , hence b/a = 2 σ q ′ exactly. Therefore, one should return RN(2 σ q ′ );

2. if q ′ = q * and r > 0 (which implies q ′ = q + ), then q ′ overestimates b * /a * . Therefore, one should return 2 σ q ′ rounded down. This is illustrated by Figure 2; 3. if q ′ = q * and r < 0 (which implies q ′ = q -), then q ′ underestimates b * /a * . Therefore, one should return 2 σ q ′ rounded up;

4. if q ′ = q * and r > 0, then q ′ underestimates b * /a * . Therefore, one should return 2 σ q ′ rounded up. This is illustrated by Figure 3; 5. if q ′ = q * and r < 0, then q ′ overestimates b * /a * . Therefore, one should return 2 σ q ′ rounded down.

When 2 σ q ′ is not a midpoint, one should of course return RN(2 σ q ′ ).

q ′ q - q * q + b * /a * Figure 2:
The number q ′ is equal to q + . In this case, the "residual" r was positive, and-since q -< b * /a * < q + -, q ′ is an overestimation of b * /a * .

q ′ q - q * q + b * /a * Figure 3: The number q ′ is equal to q * . In this case, the "residual" r was positive, and q ′ is an underestimation of b * /a * .

Proceedings of the 16th

  Symposium on Computer Arithmetic, pages 79-86. IEEE Computer Society Press, Los Alamitos, CA, 2003.


Sophie Laurens
email: sophie.laurens@math.u-bordeaux1.fr

A general family of Perfectly Matched Layers for non necessarily convex domains

Keywords: Perfectly Matched Layers, Pseudo-riemannian manifolds, Waves propagation equations, Fictive absorbing media, Friedrichs systems AMS Subject Classification 35L05, 35L40, 35Q60, 35Q61, 65N12, 78A40, 78A45

Since introduced by Berenger [1], the Perfectly Matched Layers method (PML) has become a popular approach for non reflecting Absorbing Boundary Conditions (ABC) in the numerical solution of the Maxwell equations on unbounded domains. However, most of the formulations only concern parallelepiped or simply convex domains. The goal of this paper is to present the theory for non necessarily convex PML on systems like Maxwell's.

Introduction

Whenever one solves a Partial Differential Equations numerically (by a volume discretization), one has to truncate the domain in some way. A key question is how to terminate the mesh without creating excessive echoes from the artificial truncation surface that may spoil the quality of the solution. In 1994, Berenger changed the question : instead of finding an absorbing boundary condition, he found an absorbing boundary layer. That is an artificial material independent of the boundary condition. When a wave enters the absorbing layer, it is attenuated by absorption and decays exponentially. Even if it reflects off the boundary, the returning wave after one round trip through the absorbing layer is exponentially small. Moreover, waves do not reflect at the interface. Although PML were originally derived from electromagnetism and Maxwell equations, the same ideas can immediately be applied to other waves equations.

Chew and Wheedon [START_REF] Chew | A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF] introduced the notion of complex coordinates stretching, based on analytic continuation of Maxwell's equations into spatial complex coordinates where the fields are exponentially decaying. In this paper, the stretching is interpreted as writing the same equations in a complex tangent bundle of a flat manifold in C n . For a general system, two cases appear. The first one concerns equations with an intrinsic form given by exterior derivatives and Hodge operators. In that case, we can directly study their integral inverses and thus obtain a general PML theory. This allows to reduce the non necessarily convex PML domain to the closest of the real field region. This leads to a dispersive formulation : the spatial operator is local and the time one is pseudo-differential but a localization can be performed thanks to extra Ordinary Differential Equations. On the contrary, in a second case, like for example aeroacoustics' equations with convection, such a formulation is impossible to obtain. This brings difficulties to prove the existence and uniqueness of the solution, but it is possible to prove that an exponentially absorbing and unreflected solution exists in the entire domain.

This paper is about to generalize the PML formulation for non necessarily convex domains throught a diffeomorphism ϕ, defined thanks to an affinity normal to the convexified boundary, and a strictly incresing function f , whose properties will be discussed. Such a formulation will also included "classical" PML, such as cartesian and convex ones.

Flat complex manifolds

In 2001, Lassas et al. [START_REF] Lassas | Complex Riemannian metric and absorbing boundary conditions[END_REF] showed that all PML can be obtained through a complexification of coordinates that corresponds to flat complex manifolds. Instead of stretching the coordinates, they changed the metric defined on R 3 . This method presents several advantages. First, when Maxwell's equations are written in terms of 1-forms, the differential operators take form of exterior derivatives. Second, the stretching of the metric allows to treat more general scattering geometries than before. Finally, this formulation is completely invariant as it is done without a reference to specific coordinate systems.

The manifolds used in this paper are called pseudo-riemannian in C 3 , meaning real of dimension 3, with complex tangent and cotangent bundles, with a symmetric complex metric g(., .) such as its determinant g = det(g(., .)) can always be defined with the same determination. The purpose of this section is to geometrically describe the complexification of coordinates that defines the PML.

Let (M, g) be a real manifold with a complex tangent bundle T M and a riemannian metric g jl . The matrix G = (g ij ) is associated to the metric g(., .). The tangent bundle set T x M for all x in M, is the sum U + iV , with U and V real vectors tangent to x and its complex dual is the cotangent bundle M . An element of T x M (a vector) is a complex linear combinaison on this basis

j X j (x) ∂ ∂x j . The elements of T * x M are covectors on the form j X j (x)dx j , with X j (x) ∈ C. The bilinear application T M × T M → T M : (u, v) → ∇ u v with the properties ∇ u (v + w) = ∇ u v + ∇ u w ∇ f u v = f ∇ u v and ∇ u (f v) = u(f ) v+f ∇ u v is a linear connection.
In particular, the Levi Civita connection is the only linear connection that conserves the metric g(., .) and can be expressed through the Lie brackets

[u, v] of u and v. A connection is torsion free if ∇ u v -∇ v u-[u, v] = 0 and ∇ i g = 0. The complex curvature is R(u, v)w = ∇ u ∇ v w -∇ v ∇ u w -∇ [u,v]
w. On a basis of T M , as ∂ x i , ∂ x j = 0, the nullity of the curvature means the commutation of the covariant derivatives. If η and τ are two p-forms, and if η¬τ = η J τ J is their scalar product (J stands for the set of the ordonate indices p), the Hodge star operator ⋆ is given by the relationship :

η ∧ ⋆τ = τ ∧ ⋆η = (η¬τ ) (g) dx 1 ∧ • • • ∧ dx n .
Moreover, if n is the dimension of the manifold, then ⋆ -1 τ = (-1) p(n-p) τ . The canonic euclidian complex metric g C (., .) in C 3 will be considered as the real submanifold of dimension 3 of R 6 with complex tangent and cotangent bundles, and used in cartesian coordinates

∀(u, v) ∈ C 3 , g C (u, v) = 3 j=1 u j v j . If x → x(x) is a change of variables, a complex metric on R 3 , denoted g x : T x M R 3 × T x M R 3 → C 3 , can be defined by g x (u, v) = g C (d x(u), d x(v))
, where d x is the differential of x.

In the natural frame of reference, the tangent vectors always commute but the covariant derivates only commute if the manifold is flat, which provides an intrinsic way of writing the determinant and the symbol of the operator. A flat manifold therefore defines a parallelism structure that generalizes the notion of convolution.

Formalism

The spatial domains involved are : a bounded smooth obstacle O, non necessarily convex that can be empty, a connex domain Ω, non necessarily simply connex that represents the domain of interests, δΩ 0 is the boundary between Ω and O, δΩ ∞ is the external boundary. Both are C 1 , and δΩ ∞ can be rejected to infinity if Ω is not bounded.

Consider the real manifolds with values in C 3 described by the change of variables, with p = ε + iω ∈ C,

x = ϕ(x) + 1 p f (ϕ(x)). ( 1
)
Hypotheses 1. The couple (ϕ, f ) belongs to C 1 (R 3 , R 3 ). ∀x ∈ Ω, we have ϕ(x) = x and f (x) = 0. Moreover, |ϕ(x)| = O |x|→∞ (|x|)
, f is strictly increasing with a linear growth at infinity and it exists S a convex function such as f = gradS ϕ .

These hypotheses will be justified later.

Remark : the regularity C 1 pcw (R 3 , R 3 ) ∩ C 0 (R 3 , R 3
) for (ϕ, f ) could be enough if the edges are lipschitzian (no jump discontinuities) and it could be possible to write the weak exterior differentials thanks to the unit partition theorem.

The interest of such a formulation is that the complexification of each component is made through a single function given on all R 3 and does not dependent anymore of the coordinates. Moreover, it does not take into account the shape of the PML or their absorbing direction. As f | Ω = 0 Ω and ϕ | Ω = id Ω , the system of equations is unchanged inside the studied domain. If the functions are C k , then the manifold M defined by ( 1) is a submanifold of class C k and dimension 3 in the R vectorial space C 3 of dimension 6.

Remark : as announced, this formulation is very general : chosing ϕ = x and f (x) = (f i (x)) with

f i (x) =
x 0 σ i (τ )dτ and σ i a positive function that tends to a constant at infinity, means a cartesian formulation, while f i (x) = ∂ i h(dist(x, ∂Ω)), with h an increasing convex function asymptotically linear is a convex formulation. In both cases, f is the gradient of a convex function and respects Hyp.

(1)

Helmholtz problem

Let (M, g) be an pseudo-riemannian manifold. We define the Hodge star operator ⋆ corresponding to the complex metric g, with U ∧ ⋆V = g(U, V )dvol g and dvol g = d x 1 ∧ d x 2 ∧ d x 3 . The functions x j are the components of the embedding ∼ : M → C 3 that defines the immersion. The generalized Laplacian ∆ r g for r-forms on the manifold M can be generalized for to the metric g(., .) through ∆ r g = (-1) r (⋆d ⋆ dd ⋆ d⋆). If the manifold is flat, [∇ i , ∇ j ] = 0, the components of 1-forms are given by the metric and ∆ 0 g corresponds to the opposite of the Laplace-Beltrami operator ∆ : ϕ → g ij ∇ i ∇ j ϕ.

Definition 1. (M, g) is an absorbing pseudo-riemannian manifold if 1. the manifold (M, g) is flat and M is diffeomorphic to R 3 (with φ the diffeomorphism), 2. it exists Ω ⊂ M a relatively compact open set where the metric is real and euclidian and given by g = φ ⋆ g e , with φ ⋆ the pull-back of φ and g e the euclidian metric,

3. ∀v ∈ T R x , v = 0, g x (v, v) = 0,
4. the immersion ∼ : M → C 3 guarantees that the imaginary part of x(x 1 )x(x 2 ) is signed and of constant sign for all (x 1 , x 2 ) ∈ M 2 .

For now on, the manifold considered (M, g) will be an absorbing pseudo-riemannian manifold. The fourth point of Def. (1) uniquely determines the square root

{ x -y} =   3 j=1 ( x j -y j ) 2   1/2 .
If { x -y} ≥ 0 is positive, M is called an outgoing absorbing pseudo-riemannian manifold.

Theorem 1. The fundamental solution G H for the g-Helmholtz operator for 0-forms on (M, g)

∀y ∈ M, (∆ g -p 2 )Φ(., y) = -δ y is G H (x, y) = exp(-i ω { x -y}) 4π { x -y}
with G H the usual solution for the euclidian metric, δ y a 0-current. When |x -y| → ∞, G H and ∇G H have an exponential decay in exp (-O(|x -y|)) (with p = iω + 0)

Remark : the Dirac delta is to be interpreted with respect to the volume form defined by the

metric g : if ψ is a C ∞ 0-form on M, M ψ(x)δ y (x)dvol g (x) = ψ(y).
The proof of this theorem can be found in Lassas [START_REF] Lassas | Complex Riemannian metric and absorbing boundary conditions[END_REF]. For the fundamental solution of the Helmholtz equation to be defined, the fourth condition of Def. (1) must be satisfied. The change of variables is given by Eq. (1), and by considering the limit case p = iω, we have

{ x -y} 2 = [ϕ(x) -ϕ(y)] 2 - 1 ω 2 [f (ϕ(x)) -f (ϕ(y))] 2 - 2i ω [ϕ(x) -ϕ(y)] . [f (ϕ(x)) -f (ϕ(y))] . Let us denote [f ] = [f (ϕ(x)) -f (ϕ(y))] and [ϕ] = [ϕ(x) -ϕ(y)].
The real and imaginary parts must not be null simultaneously. The condition [f ] = 0 is straight forward deduced from the previous equation. But

[f ] = 1 0 (1 -t) Df Dϕ ϕ 1 + t (ϕ 2 -ϕ 1 ) dt. If M(ϕ 1 , ϕ 2 ) is the previous matrix, a condition for [f ] = 0 is (M [ϕ] , [ϕ]) = 0, meaning that [ϕ]
is a eigenvector of the matrix M + M T associated to a null eigenvalue. A sufficient condition is M positive or null. For a continuous determination of the square root, the imaginary part must be of constant sign. A sufficient condition is M symmetric. In that case, the matrix Df Dϕ is also symmetric and there exists a convex function S such as f = grad S ϕ . This justifies Hyp. [START_REF] Berenger | A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF].

The study of these equations will first be made on the whole manifold M to guarantee the existence and uniqueness of the solution. Then, the problem on unbounded domain will be restrained to a bounded one with adapted boundary conditions.

If the hypotheses from the Rauch's theorem [START_REF] Rauch | Symmetric positive systems with boundary characteristic of constant multiplicity[END_REF] are satisfied, every strong solution is a weak solution, and strong solutions are smooth with second members of the same regularity. If the system has a smooth strong solution, then it has a weak solution. As a consequence, in this paper, the electric and magnetic fields e and h are exclusively smooth 1-forms. Given a metric g, there is a well-known one-to-one correspondence between vector fields and 1-forms. Let j and m be 2-forms standing for the second members (or Right Hand Side (RHS) of the equation) at coefficients in C ∞ with a compact support. The harmonic Maxwell's equation in vacuum are

p ⋆ e -d h = -j p ⋆ h + d e = -m, (2)
with p = ε + ik ∈ C, k = ω/c = ω √ ε 0 µ 0 ,
M = 0 d -d 0 = -M T ,
and ⋆(e, h) is the Hodge transformation of (e, h) defined by (⋆e, ⋆h). The Maxwell equations on intrinsic form are

(p ⋆ +M) e h = -j -m . ( 3
)
Lemma 1. If ϕ is a 0-form and A a 1-form with ∇A = 0 then ∆ 1 (ϕA) = (∆ 0 ϕ)A.

Proof. If A is a 1-form then A is a linear combinaison of d x i whose Hodge transformations are

⋆d x 1 = d x 2 ∧ d x 3 ⋆ d x 2 = d x 3 ∧ d x 1 ⋆ d x 3 = d x 1 ∧ d x 2 .
Without a loss of generality, we assume that A = d x 1 and

∇ i : = ∇ ∂ x i , then ∆ 1 (ϕA) = -(⋆d ⋆ d -d ⋆ d⋆)(ϕd x 1 ) = (-∂ 11 ϕ -∂ 22 ϕ -∂ 33 ϕ)d x 1 + (∂ 12 ϕ -∂ 12 ϕ)d x 2 + (∂ 13 ϕ -∂ 13 ϕ)d x 3 = -(∆ 0 ϕ)d x 1 .
Similar results are obtained with A = d x 2 and A = d x 3 .

Theorem 2. If A and B are 2-forms such as ∇A = ∇B = 0, the application G A,B (x, y) (defined thanks to the Green function G H (x, y) provided by Thm. 1)

G A,B (x, y) = p -p -1 d ⋆ d⋆ -⋆ d ⋆d p -p -1 d ⋆ d⋆ G H (x, y)A G H (x, y)B .
satisfies the properties of the Green function for Maxwell's equations on (M, g) and

∀y ∈ M, p + ⋆M T G A,B (., y) = Aδ y Bδ y . ( 4
) Proof. As p + ⋆M T = p ⋆d -⋆ d p , it comes p + ⋆M T G A,B (., y) = p ⋆d -⋆ d p p -p -1 d ⋆ d⋆ -⋆ d ⋆d p -p -1 d ⋆ d⋆ G H (., y)A G H (., y)B = p 2 -d ⋆ d ⋆ + ⋆ d ⋆ d 0 0 p 2 -d ⋆ d ⋆ + ⋆ d ⋆ d G H (., y)A G H (., y)B = (p 2 -∆ g )G H (., y)A (p 2 -∆ g )G H (., y)B = Aδ y Bδ y (lemma (1))
3.1 Case of the Maxwell equations with RHS = 0 and O = ∅

For the rest of this paper, we assume that solutions are smooth enough without loss of generality thanks to the equivalence provided by the Rauch Theorem [START_REF] Rauch | Symmetric positive systems with boundary characteristic of constant multiplicity[END_REF] between strong and weak solutions and the Friedrichs systems finally obtained. The currents will be written as differential forms with distributional coefficients since equations are intrinsic.

Theorem 3. If (M, g) is an outgoing absorbing pseudo-riemannian manifold, then the problem given by (2) without any scattering object has a unique solution

e i = M G d x i ,0 ∧ -j -m h i = M G 0,d x i ∧ -j -m . Proof. With X, Y 1-forms, we denote ⋆ X Y = ⋆X ⋆Y . Their exterior product, with P, Q 1-forms, is X Y ∧ P Q = X ∧ P + Y ∧ Q. Moreover, P * ∧ Q = -P ∧ Q * .
The weak formulation of Maxwell's equations, with A and B 1-forms, is

(p ⋆ +M) e h ∧ ϕ(A, B) = e h ∧ (-p ⋆ ϕ(A, B) + Mϕ(A, B))
By applying ⋆ to Eq. ( 4), we have p ⋆ +M T G A,B (., y) = ⋆Aδ y ⋆Bδ y . Then,

e ∧ ⋆A + h ∧ ⋆B = M e h ∧ p ⋆ +M T G A,B (., y) dvol g (y) = M e h ∧ (p ⋆ G A,B (., y)) dvol g (y) + M e h ∧ M T G A,B (., y) dvol g (y) By commutation of ⋆, M e h ∧ (p ⋆ G A,B (., y)) = M G A,B (., y) ∧ p ⋆ e h . If α and β 1-forms, d (α ∧ β) = d α ∧ β -α ∧ d β, therefore 0 = M d (α ∧ β) = M d α ∧ β - M α ∧ d β,
and then

M d α ∧ β = M α ∧ d β.
With the properties of the Green function,

M e h ∧ M T G A,B (., y) = M M T e h ∧ G A,B (., y) = M G A,B (., y) ∧ M e h .
Finally

e ∧ ⋆A + h ∧ ⋆B = M G A,B (., y) ∧ (p ⋆ +M) e h dvol g (y) = Eq.(3) M G A,B (., y) ∧ -j -m dvol g (y).
This formula is established for any 1-forms A and B. As e = e i d x i and h = h i d x i , by choosing respectively A = d x i , B = 0 and A = 0, B = d x i , the components e i and h i of the electric and magnetic fields are determined.

e ∧ ⋆A + h ∧ ⋆B = ∂Ω∞ e ∧ G 2 -h ∧ G 1 dvol g (y).
Proof. With A and B 1-forms,

e ∧ ⋆A + h ∧ ⋆B = R 3 /0 e h ∧ p ⋆ +M T G A,B (., y) dvol g (y) = R 3 /0 e h ∧ (p ⋆ G A,B (., y)) dvol g (y) + R 3 /0 e h ∧ M T G A,B (., y) dvol g (y)
Eq. ( 3) with RHS= 0 becomes

p ⋆ e h = -M e h . But -M = M T , so R 3 /0 e h ∧ (p ⋆ G A,B (., y)) = R 3 /0 G A,B (., y) ∧ p ⋆ e h = R 3 /0 G A,B (., y) ∧ 0 d -d 0 e h .
Moreover,

R 3 /0 e h ∧ M T G A,B (., y) = R 3 /0 e h ∧ 0 d -d 0 G A,B (., y) . With G A,B (., y) = G 1 G 2 , we have e ∧ ⋆A + h ∧ ⋆B = R 3 /0 e ∧ dG 2 + G 2 ∧ de -h ∧ dG 1 + G 1 ∧ (-dh) dvol g (y) = R 3 /0 e ∧ dG 2 -de ∧ G 2 -h ∧ dG 1 + dh ∧ G 1 dvol g (y) = R 3 /0 d(e ∧ G 2 ) -d(h ∧ G 1 ) dvol g (y). With Stokes formula, e ∧ ⋆A + h ∧ ⋆B = ∂Ω∞ e ∧ G 2 -h ∧ G 1 dvol g (y)
.

Remark : this theorem gives the solution for a scattering problem for the Maxwell equations, but also an estimation of the error committed by bounding the domain. In all the space, the restriction of the solution coincides with the real solution. If the problem is well-posed, there is a correspondance between the inhomogeous problem and the trace of the solution. If the solution is homogeous, this formula gives the error. Moreover, if the errors created on the articifial boundary are small, it is useless to set an artificial boundary far away for the studied domain D because of the exponential decay of the Green function : the traces of the PML solution are small and exponentially decay while they return inside D. Finally, this formula allows the superposition of solutions : for a perfect scatterer, the solution of Maxwell's equations depends on the trace of the fields on this object's boundary, which leads to the following corollary.

Corollary 5. On (M, g), Problem (2) has a unique solution

E = ∂Ω 0 -η ∧ G d x j ,0 2 + γ ∧ G d x j ,0 1 d x j , H = ∂Ω 0 -η ∧ G 0,d x j 2 + γ ∧ G 0,d x j 1 d x j .
with η and γ smooth enough currents.

Harmonic Problem

In this section, we prove that the PML problem for Maxwell's and waves equations can be written on the form

K(p, x) e h + 0 -∇× ∇× 0 e h = -j -m
with (j, m) smooth enough, as well as the following theorem.

Theorem 6. With the following hypotheses :

• K(p, x) L ∞ in x • K(p, x) holomorphic in p ∈ H ⊃ iR + * • ∃p 0 ∈ H such as Re(K(p 0 , x)) is coercive • ∃p 1 ∈ H such as K(p 1 , x) + A c i ∂ i with a compact resolvent,
for all p ∈ H/S where S is a locally finite set, eventually empty, and

Ω ∞ bounded, (K(p, x) + A c i ∂ i ) -1 is bounded in L 2
(Ω) m with maximal monotone boundary conditions.

Intrinsic form of equations

The contravariant coordinates of a random point x are denoted by (x 1 , x 2 , x 3 ), the basis vector of the cotangent bundle by dx i and dx i stands for the elementary 2-form given by the exterior product dx j ∧dx k with j = i, k = i and j < k. With E and H 1-forms, (M, g) an absorbing pseudo-riemannian manifold of metric g(., .), the Maxwell equations are the following equalities of 2-forms

p ⋆ E -dH = 0, p ⋆ H + dE = 0.
With G matrix of the metric and g = det(G), by the Hodge star operator ⋆ on M, a 1-form E = e i dx i is changed onto a 2-form whose component in dx i is √ gg ij e j (g ij components of the inverse of g(., ), √ g complex square root of the determinant g). Therefore ⋆E = √ gG -1 e. The intrinsic form of Maxwell equations is

p √ gG -1 E -dH = 0, p √ gG -1 H + dE = 0.
The purpose of this paper is to restrain the solution of Maxwell's equations to a bounded domain Ω of R 3 . The first step is to write the equations in a more appropriated system of coordinates. Let us remind that the PML media is defined by a stretching of coordinates decribed by [START_REF] Berenger | A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF]. If J(p, x) is the jacobian matrix of this transformation, then

G(p, x) = J(p, x) T J(p, x).
The associated Hodge transformation is

⋆ = √ gG -1 (p, x) = K(p, x) = det(J(p, x))J(p, x) -1 J(p, x) T -1 . ( 5
)
The domain of the family of operator (K(p, x) + A) p∈D 0 is holomorphic and independent of the frequency ω. As solving a PDE means to inverse its operator, it seems natural to look precisely the kind of operator we have to deal with. The book of Kato [START_REF] Kato | Perturbation theory for linear operators[END_REF] gives various tools to prove existence and uniqueness of solution as long as the family of operator is of type A, which is the case here.

Main results on operators of type A

Definition 2. Let X,Y be two Banach spaces. A family T (u) ∈ C(X, Y) defined for u in a domain D 0 of the complex plane is said to be holomorphic of type A if D(T (u)) = D is independent of u, and

∀(u, v) ∈ D 0 × D, T (u)v is holomorphic.
Theorem 7. Let T (z) ∈ C(X) be a holomorphic family of type A on a domain D 0 of C. We suppose that the resolvent set of T (z) is not empty ∀z ∈ D 0 and that there exists a point z 0 in D 0 such as the resolvent of T (z 0 ) is compact. Then, for all z in D 0 , the resolvant of T (z) is compact.

Harmonic Maxwell equations in the PML

In the previous section, we have established that harmonic Maxwell equations in the PML media can be written as (K(p, x)

+ A(x))u(x) = f (x), with K : (p, x) ∈ D 0 × Ω -→ Hom(C 6 ) a holomor- phic function in p = ε + iω ∈ D 0 ⊂ C, in L ∞ (Ω) for x, and A = 0 -∇× ∇× 0 . We suppose that D 0 is a domain of C including a segment like [iω min , iω max ], f ∈ L 2 (Ω) and A is Maxwell's operator whose domain D(A) is independent of p ∈ C.
Solving this equation means inversing the operator (K(p, .) + A) of domain D(A) at p = iω, which is an operator of type A. Moreover, (A, D(A)) is maximal monotone and its resolvent set is not empty. As K(p, .) is L ∞ (Ω), the operator (K(p, .) + A, D(A)) has a non empty resolvent set for all p in D 0 .

Thm. [START_REF] Rauch | Symmetric positive systems with boundary characteristic of constant multiplicity[END_REF], the main point of this section, will prove that studying the (K(p, .) + A is an operator of type A with a compact resolvent set, and the study of this set's compactness can be reduced to the study in a single point, and therefore ensures the existence and unicity of the solution.

Theorem 8 (Petkov [START_REF] Petkov | Scattering theory for hyperbolic operators[END_REF]).

If Q = div 0 0 div , for all u in D(A) ∩ (H(div, Ω) 2 ), then u (H 1 (Ω)) 6 ≤ C u (L 2 (Ω)) 6 + Au (L 2 (Ω)) 6 + Qu (L 2 (Ω)) 6 . Definition 3. If P 0 : L 2 (Ω) 3 → H(div 0 , Ω) and P ⊥ : L 2 (Ω) 3 → (grad(H 1 0 (Ω)) 2
) are the projections associated to the Hodge decomposition, for all g ∈ L 2 (Ω) 3 , g 0 = P 0 g, g ⊥ = P ⊥ g. The space H(div 0 , Ω) 2 is stable under the action of A.

Lemma 2. The operator P ⊥ K(p, .)P ⊥ is invertible on (grad(H 1 0 (Ω))) 2 for all p in D 0 except for a finite subset (eventually empty) of D 0 , denoted by S. Moreover, its inverse is holomorphic on D 0 \S .

Proof. Let h be in H -1 (Ω) 2 . Inversing P ⊥ K(p, .)P ⊥ on (grad(H 1 0 (Ω))) 2 requiere to find ω in

H 1 0 (Ω) 2 such as ∀v ∈ H 1 0 (Ω) 2 , Ω K(p, .)∇ω, ∇v dx =< h, v > H -1 (Ω) 2 ×H 1 0 (Ω) 2 , (6)
with ∇ω = ∇ω 1 ∇ω 2 . If a(ω, v) = Ω K(p, .)∇ω, ∇v dx, L(v) =< h, v > H -1 (Ω) 2 ×H 1 0 (Ω) 2
, the applications a(., .) and L(.) are respectively bilinear and linear. A Cauchy-Schwartz inequality proves that a(., .) is continuous on H 1 0 (Ω) 2 × H 1 0 (Ω) 2 and L(.) is continuous on H 1 0 (Ω) 2 . The previous problem is in fact, to find ω in

H 1 0 (Ω) 2 such as ∀v ∈ H 1 0 (Ω) 2 , a(ω, v) = L(v). (7)
The Lax-Milgram theorem for p = p 0 guarantees that this equation is well-posed. The coercivity of K(p 0 , .) and the continuity of L(.) on H 1 0 (Ω) 2 leads to the inequality satisfied by the solution β ω H 1 0 (Ω) 2 ≤ h H -1 (Ω) 2 , with β > 0. Therefore, the resolvent of P ⊥ K(p 0 , .)P ⊥ is compact because of the compactness of the injection H 1 0 (Ω) in H -1 (Ω). By Theorem [START_REF] Rauch | Symmetric positive systems with boundary characteristic of constant multiplicity[END_REF], the resolvent of P ⊥ K(p, .)P ⊥ is compact for all p in D 0 . As a consequence, this operator is either singular on D 0 either invertible for all p in D 0 except for a finite subset (eventually empty) of D 0 denoted S. As Problem ( 7) is well-posed for p 0 , this operator is invertible ∀p in D 0 \S.

Let us show that operator (P ⊥ K(p, .)P ⊥ ) -1 is holomorphic on D 0 \S. By application of the closed graph theorem to the closed set grad(H 1 0 (Ω)) 2 of L 2 (Ω) 6 , we have (P ⊥ K(p, .)P ⊥ ) -1 in L((grad(H 1 0 (Ω))) 2 ). As P ⊥ K(p, .)P ⊥ is holomorphic on D 0 , its inverse is holomorphic on D 0 \S.

Theorem 9 (Kato [START_REF] Kato | Perturbation theory for linear operators[END_REF]). If it exists p 0 ∈ D 0 such as K(p 0 , .) coercive, then harmonic Maxwell system in the PML media is well-posed for all real frequency except for a finite discrete subset (eventually empty) of R.

Proof. By projection of Eq. ( 6) with the Hodge decomposition, the problem is to find u = u 0 + u ⊥ ∈ D(A) such as

P 0 K(p, .) u + A u 0 = f 0 , P ⊥ K(p, .)u = f ⊥ .
As u = P 0 u + P ⊥ u, the previous system becomes P 0 K(p, .)P 0 u + P 0 K(p, .)P ⊥ u + Au = f 0 , P ⊥ K(p, .)P ⊥ u + P ⊥ K(p, .)P 0 u = f ⊥ .

With Lem. (2), if A is the operator defined by the restriction of the operator A to the set (H(div 0 , Ω)) 2 , then D( A) = D(A) ∩ (H(div 0 , Ω)) 2 , so the problem is to inverse in (H(div 0 , Ω)) 2 the closed operator R(p) = P 0 K 0 (p, .)P 0 -P 0 K 0 (p, .)P ⊥ (P ⊥ K 0 (p, .)P ⊥ ) -1 P ⊥ K 0 (p, .)P 0 + P 0 K 1 (p, .) A

The operator R(p) is holomorphic on D 0 \S = D 1 . By use of Thm. ( 8) and ( 7) to the family of holomorphic operator R(p) + A of type A, it appears that R(p) + A is a holomorphic family of closed operators with compact resolvents for p ∈ D 1 . As a consequence, this operator is either singular on D 1 either invertible for all p in D 1 except for a locally finite discrete subset (eventually empty) of D 1 . Operator P ⊥ K(p 0 , .)P ⊥ defines a sesquilinear coercive form in (grad(H 1 0 (Ω))) 2 because K(p 0 ) is coercive in L 2 (Ω) 6 which implies P ⊥ K(p 0 , .)P ⊥ invertible. A is maximal monotone and K(p 0 , .) is monotone coercive and bounded in L 2 (Ω) 6 , thus the operator A + K(p 0 , .) is invertible. Finally, the operator R(p 0 , .) + A is invertible on (H(div 0 , Ω)) 2 . So R(p) + A is invertible for all p in D 1 except for a locally finite discrete subset (eventually empty) of D 1 . This proves that the operator K(p, .) + A is invertible on L 2 (Ω) 6 for all p in D 1 except for a locally finite discrete subset (eventually empty) of D 1 .

Study of the coercivity of the PML matrix

According to Thm. (9), to have existence and uniqueness of the solution, it must exists a p 0 in C such as K(p 0 , x) is coercive. If we exhibit such a point, the resolvent is not empty and there is no more alternative.

Property 10. The matrix K(1, x) is symmetric definite and positive for all x in R 3 .

Proof. As K(p, x) is defined by Eq. ( 5), we have

K(p, x) T = det(J(p, x))J(p, x) -1 J(p, x) T -1 T = det(J(p, x))J(p, x) -1 J(p, x) T -1 = K(p, x), with J(p, x) = D( x) D(ϕ(x)) D(ϕ(x)) D(x) . So K(p, x) is symmetric for all p ∈ C, and a fortiori K(1, x) is symmetric. If X ∈ R 3 with X = 0, then X T K(1, x)X = X T det(J(1, x))J(1, x) -1 J(1, x) T -1 X = det(J(1, x))X T J(1, x) -1 J(1, x) T -1 X = det(J(1, x))X T J(1, x) -1 X T J(1, x) -1 T . If A = X T J(1, x) -1 , then X T K(1, x)X = det(J(1, x))AA T .
As AA T is definite positive for all matrix A, we deduced that K(1, x) is definite positive if det(J(1, x)) > 0. Let F(x) = D(ϕ(x)) D(x) be the jacobian matrix of ϕ. As ϕ is an increasing function that represents the embedding in a manifold M of C 3 , we have det(F(x)) > 0. Let us denote F 0 (x) =

D( x) D(ϕ(x))

. For all (i, j) ∈ 1, n ,

F 0 (x) ij = δ ij + f ′ (ϕ(x)
). As f is strictly increasing, det(I 3 + F 0 (x)) > 0.

Explicit formulations for Maxwell's Equations

Lemma 3. The matrix of the harmonic PML formulation is

K(p, x) = detT(x) T(x) -1 P * MPT(x) * -1 with T(x) = Dϕ Dx , M = diag    i ∨ k (p + λ k )    and P * = P -1 .
Proof. Given the embedding x = ϕ(x)

+ 1 p f (ϕ(x)), it
Ī(x) ∂ ∂t e h + Ā(x) e h + B(x) u v + 0 -∇× ∇× 0 e h = -j -m ,
where (u, v) are solutions of the ODE

T(x) * P -1 0 0 T(x) * P -1 ∂ ∂t u v + T(x) * P -1 F(x) 0 0 T(x) * P -1 F(x) u v = e h , with F(x) = diag (λ i ) , A(x) = diag    i ∨ k λ k -λ i    , B(x) = diag    i ∨ k (λ i -λ k )    and Ī(x) = M(x) M(x) * -1 , Ā(x) = M(x) A(x) 0 0 A(x) M(x) * -1 B(x) = M(x) B(x) 0 0 B(x) M(x) * -1 with M(x) = det(T(x)) T(x) -1 P * 0 0 T(x) -1 P * .
Proof. Let (i, j, k) be three different indices.

(p + λ i )(p + λ j ) p + λ k = p + (λ i + λ j -λ k ) - (λ i -λ k )(λ j -λ k ) p + λ k ,
The diagonal matrix of Lem. (3) can be decomposed as

pI 3 + A(x) + B(x) (pI 3 + F(x)) -1 with F(x) = diag (λ i ) A(x) = diag    i ∨ k λ k -λ i    B(x) = diag    i ∨ k (λ i -λ k )    .
With T(x) = Dϕ Dx , the PML matrix can be decomposed as

det(T(x)) T(x) -1 P * pI 3 + A(x) + B(x) (pI 3 + F(x)) -1 PT(x) * -1 . ( 8
)
We notice that (pI 3 + F(x)) -1 PT(x) * -1 = T(x) * P -1 (pI 3 + F(x))

-1

. By choosing Ī(x), Ā(x), B(x) and F(x) as written above, the Maxwell equations in the PML become

p Ī(x) + Ā(x) + B(x) p I(x) + F(x) -1 e h + 0 -∇× ∇× 0 e h = -j -m . If u v = p I(x) + F(x) -1 e h
, it is a solution of an ODE that can be introduced in Maxwell's

equations on R 6 to have          p Ī(x) e h + Ā(x) e h + B(x) u v + 0 -∇× ∇× 0 e h = -j -m p I(x) u v + F(x) u v = e h .
Using the Laplace inverse transformation ends the proof.

Important remark : this theorem shows an important property of the PML. The spatial part of the operator 0 -∇× ∇× 0 is not modified in the PML media. This guarantees the possibility to choose without constrain the numerical approximation method (Finite Elements, Finites Volumes, Discontinuous Galerkin, . . . ). Moreover, a Gedney formulation can be obtained with an other simple element decomposition (for example for an FTDT approximation with leap frog schemes).

The previous formula has been established for Maxwell's equations in 3D. For bidimensional equations, the formula is deduced for Thm. (11) through an Hadamard's method of descent.

6 Some numerical PML examples for a L-shaped geometry

A general non convex formulation

The embedding is given by θ = θ,

ρ = φ(ρ, θ) + 1 p f (φ(ρ) -R) ,
and the diffeomorphism is φ(ρ, θ) = kρ + (1k)ρ 0 (θ). The eigenvalues of the jacobian matrix are

α(ρ, θ) = f ′ (φ(ρ, θ) -R), β(ρ, θ) = f (φ(ρ, θ) -R) φ(ρ, θ)
. The PML matrix is

K(p, ρ, θ) = Q(θ)P(ρ, θ) p I(ρ, θ) + A(ρ, θ) + B(ρ, θ) pI 3 + F(ρ, θ) -1 P T (ρ, θ)Q(θ) T , (9)
with

I(ρ, θ) = k φ(ρ, θ) ρ I 3 , F(ρ, θ) = diag α(ρ, θ), β(ρ, θ), 0 , B(ρ, θ) = k φ(ρ, θ) ρ diag α(ρ, θ)(α(ρ, θ) -β(ρ, θ)), β(ρ, θ)(β(ρ, θ) -α(ρ, θ)), α(ρ, θ)β(ρ, θ) , A(ρ, θ) = k φ(ρ, θ) ρ diag β(ρ, θ) -α(ρ, θ), α(ρ, θ) -β(ρ, θ), α(ρ, θ) + β(ρ, θ) , P(ρ, θ) =       φ(ρ, θ) ρ -2 (1 -k) ρ ′ 0 (θ) ρ 0 0 k 0 0 0 k φ(ρ, θ) ρ -1      
, and Q(θ) the rotation matrix around the zaxis. In that case, the matrix behind p is no longer the identity matrix. If [ε] corresponds to this part of the matrix det(J(x)) J(x) T J(x)

-1 , meaning A 0 = [ε] 0 0 [ε]
, the CFL condition associated is defined though λ 0 the smallest eigenvalue of [ε] by dt dx ≤ λ 0 . As its value depends on the mesh, the CFL condition is the maximal value of

k 2 + 4 (1 -k) ρ ′ 0 (θ) ρ 2 + k + (1 -k) ρ 0 (θ) ρ 2 2k k + (1 -k) ρ 0 (θ) ρ - k 2 + 4 (1 -k) ρ ′ 0 (θ) ρ 2 + k + (1 -k) ρ 0 (θ) ρ 2 2 -1 4k 2 k + (1 -k) ρ 0 (θ) ρ 2 .
A numerical visualisation on an unstructured mesh of the previous formula shows that the CFL condition can reach 180, but only for few elements, where the values of ρ ′ 0 (θ) are important. As the matrix is well-conditioned, this formulation can be used for harmonic problems. The domain is reduced by comparison with a convex formulation, and few more complexity is involved. But for unsteady problems, the penalization induced by the CFL condition may no be balanced by the reduction of the triangles' number. A first solution is to work implicitely.

A more specific formulation

A part of a L-shaped geometry is convex, so an other idea is to establish a new formula that takes advantage of this convex part. The PML domain is divided in two parts : one with a cartesian formula, and the other one with a non convex formula, that matched continously the cartesian PML at the border of each domain.

B(x) = 1 κ diag (σ x (x)(σ x (x) -σ y (x)
), σ y (x)(σ y (x)σ x (x)), σ x (x)σ y (x)) and F(x) = diag σ x (x), σ y (x), 0 . The matrix T(x) is triangular and depends on the non convex part of the domain where it is calcu-

lated. If y ≥ x, T + =   1 1 -κ 0 0 κ 0 0 0 κ   , otherwise T -=   κ 0 0 1 -κ 1 0 0 0 κ   . The smallest eigenvalue of the matrix T(x)T(x) T is 1 κ + κ -1 -(1 - 1 κ ) √ 1 + κ 2 .
As κ ≥ 1, the maximal CFL condition is 1 and does not depend on the shape of the domain. Even if the L-shaped has an important length by comparison with its width, the CFL condition remains the same, contrarily to the first formulation established. This embedding has a strong analogy with the cartesian PML and therefore presents a dissymetry in the waves absorption.

An other solution, based on the convex PML, is to choose x = φ(x) + 1 p f φ(x) -RI 3 . The PML matrix becomes K(p, x) = T T (x) pI + A(x) + 1 p + 2σ(x) B(x) T(x), where σ stands for the derivative of the function f , the triangular matrix T(x) is unchanged,

I = 1 κ I 3 , A(x) = 2σ(x) κ   0 -1 0 -1 0 0 0 0 1   and B(x) = 2σ 2 (x) κ   1 1 0 1 1 0 0 0 0   .
The numerical experiments of this paper are obtained with this formulation.

Numerical Simulations

All the following simulations are made with the Maxwell Bidimensional Equations, with Silver Müller conditions at the boudary of the PML domain, and meshes are unstructured. The origin of the reference frame is the corner of the L, bewteen the domain of study and the non convex PML.

Two kind of problems will be computed, which correspond to the cases studied in Sec. 3. The first simulation is the propagation of a RHS

f (x) = 0, 0, exp (x -x c ) 2 + (y -y c ) 2 (x -x c ) 2 + (y -y c ) 2 -R 2 1 D (x) T ,
where D is the disc of center (x c , y c ) = (-1, -1) and radius R = 0.5. The pulsation is ω = 2π, the σ function is of order 2, with a coefficient that guarantees a decrease in the PML at 10 -2 , and κ = 1/3. For the comparison with cartesian PML, the pulsation chosen is ω = 3π, the order of the sigma function is 2 and the mesh is modified to use a thicker scattering object. On a visual point of view, the non convex PML and the cartesian ones lead to the same solution in the physical L-shaped domain D. To confirm this, we take into account two errors. Let ϕ N C be the non convex PML solution and ϕ C the cartesian PML one. The first error, ε 1 , corresponds to the l 2 norm of the difference on D, and the second one, ε 2 , represents the local error on each element :

ε 1 = ω k ∈D |ϕ N C (ω k ) -ϕ C (ω k )| 2 vol(ω k ) ω k |ϕ C (ω k )| 2 vol(ω k ) ε 2 (ω k ) = |ϕ N C (ω k ) -ϕ C (ω k )| 2 .
Tab. 1 gives the error ε 1 function of the number of elements per wave-length λ. The mesh are rough, with very few layers, to place ourselves in the worst case. The error ε 1 decreases slower with the refinement of the mesh for the scattering problem. The mesh of Fig. 2(o) is very rough, especially in the area used especially for the cartesian simulation. The damping starts outside the convexified of the scatterer, the waves are brutally moved away : there is a lot of dissipation. (r) ε2 with λ/16 elements per wave length with a centered" RHS (s) ε2 with λ/8 elements per wave length (case test (j))

ε 1 λ/8 λ/16 λ/32 Scattering   5.81e -2 5.65e -2 6.67e -2     1.65e -2 1.61e -2 1.76e -2     1.031e -2 9.7e -3 1.07e -2   RHS = 0   1.28e -2 5.4e -3 5e -3     5.6e -3 5.4e -3 2.3e -3     9.72e -4 9.25e -4 3e -4  
For unsteady problems, as the PML matrix contains some terms at order -1 in p, an extra equation has to be added. The Maxwell Equations on R

3 × R + are            Ī(x) ∂ϕ(x, t) ∂t + Ā(x)ϕ(x, t) + ψ(x, t) +   0 0 -∂ y 0 0 ∂ x -∂ y ∂ x 0   ϕ(x, t) = 0 ∂ψ(x, t) ∂t + 2σ(x)ψ(x, t) = B(x)ϕ(x, t).

Conclusions

This paper presents a fully proved non convex PML theory and some numerical experiments. The aim is to obtain the most general formulation possible, thanks to the families (ϕ, f ). The hypothesis of convexity is underneath, through the diffeomorphism ϕ, and absolutely necessary to prevent singularities to appear in the Green function.

The case of the L-shaped geometry can be generalized as long as it is possible to express the diffeomorphism from an affinity normal to the boundary of the convexified domain. It is always possible if the domain is not trapping. For example, a U-shaped geometry where the boundaries are strictly parallel can not give such an affinity, but if the boundaries are "opened", ϕ can be found by using the normals to the convexified.

In the non convex PML, even if the luminuous rays are twisted (consequence of the complex metric with a different wave speed), the permittivity tensor is more complex -with anisotropic permittivitiesbut the restriction of the solution to the physical domain are numerically correct. The trace of the solution on the boundary, which almost inverses the problem, will be an excellent preconditionner for integral methods. For harmonic problems, non convex PML present two advantages : the number of cells (and unknowns) is reduced, and there is no need to inverse the matrix on every elements, which is sparse. For unsteady problems, the change of the wave speed implies a drop of the CFL conditions : this saves some memories but no calculus time. As the exageration of the CFL condition depends on a very few elements, an implicit formulation can be used to solve this problem.

The comparison with cartesian PML gives excellent results, even if the choice of the numerical approximation was very poor with a Finite Volums method. The use of methods like Discontinuous Galerkin has to be considered.

In 2001, this approach of the PML theory thanks to pseudo riemannian manifolds with complex tangent and cotangent bundles was already tried by Lassas [START_REF] Lassas | Complex Riemannian metric and absorbing boundary conditions[END_REF] for harmonic problem. There are four major differences between their work and the present one. For the restriction to a bounded domain, for harmonic problem, some frequencies can be excluded (we remind that Thm. (9) specifies for all p in H except for a locally finite eventually empty set S) Lassas and Co. do not have that restriction but their theory implies to set the artificial boundary "far enough" of the domain D. In this paper, it can be placed as closed as we want. In fact, Lassas and Co. requiere an extra hypothesis to define their absorbing pseudo riemannian manifold : the asymptotic η-euclidianity of the metric, meaning they impose a strict convexity and both the metric and the connection have to be majored asymptotically by an euclidian metric. We do not need that hypothesis. Last point : we study the unsteady problems, while their work was only about harmonic formulations. Moreover, this paper explicitely gives some formulations (diffeomorphims and functions) that leads to the unsteady formulation. As the embedding is defined in a very general way (cartesian and convex PML can be described as well), and the PML matrix obtained is composed of rational fraction that can be decomposed and gives the Friedrichs unsteady system while adding an ODE. The problem of chosing the couples (ϕ, f ) is general for every non trapping domain, but has to be written. The diffeomorphism is deductible from an affinity normal to the boundary of the convexified domain.

  and the Hodge star operator ⋆ defined by the Euclidean metric to convert 1-forms to 2-forms. The fields are rescaled for symmetry, meaning e → √ ε 0 e and h → √ µ 0 h. The Maxwell operator is

3. 2

 2 Case of Maxwell's equations with RHS= 0 and O = ∅ Theorem 4. The Green function G A,B (., y) can be decomposed as G 1 G 2 , with A and B 1-forms. If (M, g) is an outgoing absorbing pseudo-riemannian manifold, we have a Stratton-Chu formula with the solution (e, h) of Problem (2)

  is possible to have a general formulation for the PML matrix. The jacobian matrix J(p, x) is D P a unitary matrix P * = P -1 such as Df Dϕ = P * I + 1 p D P with D a diagonal matrix and λ i , i = 1, 2, 3 its eigenvalues. With T(x) = Dϕ Dx , we have J(p, x) = P * I + 1 p D PT(x). The Hodge transformation associated is ⋆ = det(J(p, x))J(p, x) -1 J(p, x) T -1 . The proof is ended by writing K(p, x) = p⋆. Theorem 11. The unsteady Maxwell equations in the PML media are

Fig. 1

 1 Fig.1explains the geometry. The domain of study D, the interior L, is represented in white while the cartesian PML are with horizontal stripes. The diffeomorphim φ transforms D by adding the black area. In this specific area, whose thickness is controlled by the parameter κ, the waves do not

Figure 1 :

 1 Figure 1: The different domains for a non convex PML formulation

  (a) Real part of Ex (b) Imaginary part of Ex (c) Real part of Ey (d) Imaginary part of Ey (e) Real part of Hz (f) Imaginary part of HzAs the cartesian PML are a considered as a reference solution for the problems studied, here are comparison with the non convex PML introduced in this paper. The simulation has the same initial condition as previously, but the center is moved to the right. For the wave to propagate correctly onto the L-shaped domain, it has to go through the non convex PML. Fig. (g) shows the cartesian domain involved, and the black area is the one added in order to apply a cartesian method. Fig. (h) represents the mesh generated by the PDETOOLBOX of MatLab (like all the others used in this paper).

  (g) The different domains for a cartesian PML formulation (h) Roughest mesh for the problem with a RHS = 0 (i) Real part of Hz -cartesian PML (j) Real part of Hz -non convex PML (k) Imaginary part of Hz -cartesian PML (l) Imaginary part of Hz -non convex PMLFor the scattering problem, the study of an incident plane wave ϕ inc reflected by a L-shaped geometry Γ of length 3 and width 1 is processed, with ϕ inc (x) = (-k y , k x , 1) T exp (ω(-k.x)). k = (k x , k y ) T gives the direction of the incident wave, and it is chosen as k = conditions at the boundary of the scatterer are ∀x ∈ ∂Γ, f (x) =∂ω k ∩∂Γ -M 1 ϕ inc (x)exp (-ωk.x), with M 1 = n ⊗ n 0 -n T 0, n the outer-pointing unit normal and ω k the element considered. The PML domain has a thickness of 1/3, the function σ is of order 1. Others parameters remain the same as the previous simulation.

  (m) Imaginary part of Hz (n) Real part of Hz

  (o) Roughest mesh for the scattering problem (p) Real part of Hz -cartesian PML (q) Real part of Hz -non convex PML

Fig. (

  Fig. (r) gives the error ε 2 with a RHS = 0, whose support's center is "centered", i.e. (-1, -1) and a mesh of λ/16 elements per wave-length and Fig. (s) for a mesh of λ/8 elements per wave-length, corresponding to the test (j).

  with Ī(x) = T(x)I(x)T(x) T , Ā(x) = T(x)A(x)T(x) T and B(x) = T(x)B(x)T(x) T , with the same kind of boundary conditions as the harmonic problem :The boundary conditions at the boundary of the scatterer are ∀x ∈ ∂Γ, ∂ω k ∩∂Γ -M 1 ϕ inc (x, t)exp (-ωk.x) exp(iωt). For the simulation, an implicit Euler scheme is used. The time step is determined by ∆t = min e vol(ω e ) vol(∂ω e ) .

( t )

 t Solution of the unsteady problem : third component Hz(x, T ) with T = 1.8 Remark : the previous numerical experiments have an unnecessary domain of calculus inside the non convex PML. At the external boundary of the cartesian PML, the solution is vanishing. This property is verified in the non convex PML on a parallel of the exhaustion. As a consequence, for a L-shaped geometry, the non convex area required is less than a L. On Fig. (u) and (v), the vertical stripes represent the useless domain : the non convex PML domain is very thin. Fig. (u) corresponds to Fig. (1) while Fig. (v) has an elarged PML domain to show with more visibility the effective non convex PML domain, which is very thin by comparison with the thickness of the cartesian area required on this specific geometry.

Table 1 :

 1

ε 1 for (E x , E y , H z ) depending of the number of elements per wave length


S Topoliantz

J.-F Ponge

& P Lavelle

Humus components and biogenic structures under tropical slash-andburn agriculture

Slash-and-burn cultivation in the humid tropics can cause changes in the composition of topsoil, depending on the duration of the fallow. We studied differences between practices, using the small-volume micromorphological method, to quantify the distribution of solid components in the topsoil, concentrating on plant organs and biogenic structures created by soil animals.

We compared samples of topsoil from five plots, two at Maripasoula, an Aluku village along the river Maroni (French Guiana), with short fallow (≤8 years), and the other three at Elahe, a Wayana village along the same river, with long fallow (≥ 25 years). At both sites structures created by arthropods other than ants gave way to ones formed by ants and annelids under the influence of fire and cultivation. This change was more abrupt under long fallow, because of the time needed to restore the arthropod community. Charcoal and charred plant material were incorporated by earthworms into the mineral soil, forming dark grey to black aggregates. Charcoal became mixed with the mineral soil faster at Elahe than at Maripasoula, where it accumulated in the topsoil. The reason seems to be an imbalance between charcoal inputs (from repeated fires) and the capacity of burrowing animals (earthworms, ants) to mix it with the mineral soil.

Introduction

Traditional shifting agriculture in the tropics balanced food production and nature conservation. In French Guiana, this long-lasting equilibrium is threatened by increasing demographic pressure and settling of cultivators for social welfare, schools, medical care and administration. In Maripasoula (population 1200), a settlement along the Maroni river (bordering French Guiana and Suriname) of mostly Aluku people (of African lineage), the duration of fallow in the slash-and-burn system has decreased from 15 to 7-8 years in the last thirty years. It is much shorter than the 15 to more than 100 years in the traditional shifting cultivation still practised by Wayana Amerindians as in the village of Elahe [START_REF] Fleury | Les populations du Haut-Maroni et le projet du Parc National de la Guyane[END_REF].

Under traditional shifting cultivation (alternating short periods of cropping with long fallows), burning of woody vegetation before cultivation fertilizes the soil with ashes. During the fallow, the nutrient status of the soil, impoverished by crops, regenerates after several years from the organic matter added by the regrowth of vegetation, and we know from many accounts how important soil organic matter is for the maintenance of soil fertility in the tropics. The stabilization of soil organic matter within biogenic structures such as earthworm faeces help to preserve a reservoir of water, nutrients and space in which plant roots, animals and microbes can grow. Among other approaches, the fate of charcoal, from incomplete combustion of wood, has been considered of potential use for alternative systems of tropical agriculture [START_REF] Glaser | Ameliorating physical and chemical properties of highly weathered soils in the tropics: a review[END_REF]Lehmann et al., 2003;Glaser & Woods, 2004). [START_REF] Glaser | The Terra Preta phenomenon: a model for sustainable agriculture in the humid tropics[END_REF] found that charcoal is a source of stable and fertile humus in Amazonian Terra Preta. Our own experiments showed that the combined use of charcoal and manioc peel could sustain legume production on rather infertile tropical soils (Topoliantz et al., 2005).

The present study was done in southern French Guiana. Its aim was to compare the impact of slash-and-burn cultivation on humus components and biogenic structures in the soil of two agricultural systems, differing in the duration of the fallow (8 years in Maripasoula and >100 years in Elahe) and that of continuous cultivation (1 year in Maripasoula, 3 years in Elahe). Manioc (Manihot esculenta Cranz) was the major crop in both villages. These sites were chosen as typical examples of (i) traditional shifting agriculture, and (ii) change to permanent agriculture, on the basis of preliminary investigations on the impact of agricultural practices on soil fertility along the Maroni river [START_REF] Grandisson | Diagnostic de la reproductibilité de la fertilité dans les abattis de la région Ouest-Guyane: étude des stratégies de gestion et analyse des trajectoires d'évolution des composantes de la fertilité en fonction du milieu et des modes de mise en valeur[END_REF]. The crop (manioc) was the same and the soils (Oxisols) are similar. The two sites differ only in the way they are managed.

Materials and methods

Study sites

The first site, typical of traditional slash-and-burn shifting agriculture, is upstream of the Amerindian Wayana village of Elahe, on the Tampock river, a subsidiary of Maroni (N 3°26'; W 53°59'). Three years ago a family had cut and burned contiguous fields (abattis) within a secondary forest. This forest seemed, from a census of woody species typical of mature forests such as Astrocaryum sciophilum (Miquel) Pulle and Dicorynia guianensis Amsh. [START_REF] Poncy | The lowland high rainforest: structure and tree species diversity[END_REF], to have been let untouched for at least 100 years. We sampled an abattis which had been under cultivation for 3 years (EA), and the nearby untouched old secondary forest (EF) in July 1999, approximately 100 m from EA. The plot in the secondary forest was resampled in May 2000 after it was burnt in December 1999 for cultivation (EFB).

The second site is by the Maroni river (N 3°39'; W 54°2') near Maripasoula, 25 km downstream of the first site. It is in a large village inhabited mostly by Aluku people (black people of African lineage, descended from slaves who escaped at the end of the 18th century). During the last thirty years the increase of the population decreased the surface of cultivable land per family and the duration of the fallow [START_REF] Fleury | Les populations du Haut-Maroni et le projet du Parc National de la Guyane[END_REF]. We sampled a one-year-old abattis at the end of the crop period, opened by an Aluku family by burning an 8-yearold woody fallow (MA), and the nearby unburnt fallow (MF), approximately 100 m from MA. The latter plot was characterized by typical pioneer woody species such as Cecropia latiloba Miq. and Inga capitata Desv. Unlike the inhabitants of Elahe village, who can exploit the same abattis for two or three years because cultivation occurs after long fallow, cultivators living in Maripasoula do not crop the same field for more than one season, due to a rapid decrease of yield and the invasion of their fields by weeds and crop parasites [START_REF] Fleury | Les populations du Haut-Maroni et le projet du Parc National de la Guyane[END_REF]. As above we sampled both plots in July 1999. We intended to resample in 2000 the woody fallow plot (MF), which should have been burnt in December 1999 by the Aluku family, but as it was not burnt and so we could not resample it as for EFB.

The mean size of abattis is 2 ha at Maripasoula and 1 ha at Elahe. At both sites, manioc is the basic food crop, and is supplemented by fishing and hunting. Abattis are roughly circular, originating from cutting and burning the forest in places accessible by foot or by canoe. The people of both communities prefer soils that are sandy. The Wayanas also prefer dark soils. Both people avoid valleys. In each abattis all tree trunks and saplings are cut, except those kept as holy trees. The felled trees are allowed to dry. Some wood is taken for timber and cooking, and the rest is burned and its residues left on the field. Manioc cuttings, from previous crop, are planted after resprouting. The soil is not cultivated, but excavated locally to plant the cuttings and then refilled with a hoe. The aspect of the site is used to select places proper for burying manioc cuttings. Hollows in the undulating terrain are often used for bananas or rice.

The main differences between Aluku and Wayana practices lie in the duration of the crop. Each year Alukus cut and burn new abattis around the village, whereas the Wayanas use the same abattis for two or three years, the precise duration depending on the soil's fertility and the spontaneous establishment and resprouting of vegetation, then shift to another place. They weed only by hand during the first four months following manioc planting, and they do not use herbicides, pesticides or fertilizers.

The climate is warm (mean annual temperature 26°C) and rainy (2000 mm per year), with a main dry season from September to December and a shorter one from March to April. Table 1 lists the main physico-chemical properties of the topsoil (Oxisol) which was sampled in the study plots. At both sites the soil was sandy and acidic, although slightly less acidic in Maripasoula (pH5.0) than at Elahe (pH4.7). The soil was air-dried before transport to the laboratory for chemical analyses. Soil pH was measured electrometrically on a 1:2.5 soil:water mixture. Total C and N were measured by the dry combustion method after hydrochloric dissolution of carbonates according to ISO 10694 and ISO 13878, respectively (Anonymous, 1999).

Soil micromorphology: the small-volume method

We studied the distribution of humus components and biogenic structures in topsoil by small-volume micromorphology sensu [START_REF] Bernier | Humus form dynamics during the sylvogenetic cycle in a mountain spruce forest[END_REF], further refined for biogenic structures by [START_REF] Peltier | Humus forms in Mediterranean scrublands with aleppo pine[END_REF] and adapted by us to agricultural soils [START_REF] Topoliantz | Earthworm and enchytraeid activity under different arable farming systems, as exemplified by biogenic structures[END_REF]. This optical method allows one to identify and estimate the proportion of solid components in successive layers of a given soil profile. Combined with multivariate methods, it can be used to compare soil profiles along gradients [START_REF] Peltier | Humus forms in Mediterranean scrublands with aleppo pine[END_REF][START_REF] Frak | The influence of altitude on the distribution of subterranean organs and humus components in Vaccinium myrtillus carpets[END_REF][START_REF] Sadaka | Climatic effects on soil trophic networks and the resulting humus profiles in holm oak (Quercus rotundifolia) forests in the High Atlas of Morocco as revealed by correspondence analysis[END_REF] and in vegetation patchworks [START_REF] Patzel | The heterogeneity of humus components in a virgin beech forest[END_REF].

Five samples for micromorphological analyses and five samples for physico-chemical analyses were taken in each of the five study plots. They were regularly spaced along a 30-m transect crossing the centre of each plot. This allowed us to embrace within-plot variation while avoiding edge effects.

For micromorphological investigations we took soil blocks 7 cm x 7 cm x 6 cm (length x width x height) which we shaped with a sharp knife without disturbing litter nor soil structure. Then we separated the top 5 cm by hand in successive layers 0.5 to 3 cm thick according to their appearance, which we considered homogeneous. Each layer was immediately fixed in 95% ethyl alcohol then transported to the laboratory for further analysis. In the immediate vicinity of each sampling plot another sample (10 cm deep) was taken, air-dried as soon as possible before transport to the laboratory for physico-chemical analyses. A total of 251 layers (for micromorphological analyses) and 25 topsoil profiles (for physico-chemical analyses) were sampled.

The material from each layer was gently spread in a Petri dish (150 mm diameter), with as little disturbance of the aggregates as possible, and then covered with alcohol. Ethyl alcohol precipitates colloids, thereby helping to preserve aggregates, provided they are not crushed by forcing them with an instrument. We observed each layer under a dissecting microscope after covering the layer of solid matter with a transparent 600-pt grid, which we had previously prepared by piercing a transparency film at nodes of a 5-mm grid.

Following dot lines under the dissecting microscope and changing the focusing plane every time a dot was encountered, each element which was located just below a node of the grid was identified then counted. This method allowed us to estimate the relative volume of components of the soil matrix, including plant organs at varying stages of development (for subterranean organs) or decomposition, mineral particles of varying size and nature, aggregates of varying colour, size and shape. The Munsell® code for soil colours was used to classify colours of aggregates into five broad classes (light, light brown, brown, light grey, grey, black).

Data analysis

Micromorphological data were analysed by correspondence analysis, CA [START_REF] Greenacre | Theory and Applications of Correspondence Analysis[END_REF]. The matrix analysed crossed 81 layers (as columns) and 82 classes (as rows). Additional (passive) variables were added as rows, in order (i) to facilitate the interpretation of the factor axes, and (ii) to discern mean trends in the vertical distribution of humus components in the five plots studied.

They comprised the five plots (MA, MF, EA, EF, EFB) and, for each plot, five depth levels (0-1 cm, 1-2 cm, 2-3 cm, 3-4 cm, 4-5 cm).

Results

Table 2 lists the 82 classes that were identified under the dissecting microscope. The projection of classes (topsoil components) in the plane of the first two factor axes of the CA (11% and 8% of total variance, respectively) displayed three branches (Figure 1). The percentages of variance extracted by axes 1 and 2 (19%) far exceeded the confidence interval given by random permutation of rows and columns [START_REF] Lebart | Traitement des Données Statistiques. Méthodes et Programmes[END_REF] At some depth (varying from 2 cm in MA to 3 or 4 cm in other plots) all sampled layers exhibited similar features, factor coordinates being negative both for axes 1 and 2. Prominent differences between plots were exhibited mostly in the top 1 cm of the soil. The composition of the top centimetre in MA did not differ to a great extent from that of forest soils (MF and EF), both showing positive values along axis 1, although differences between surface and deeper soil were much less pronounced in MA, because there was less litter (Figure 2). In contrast, at Elahe the top 2 cm of the soil was strongly affected by the shift to agriculture (EA, EFB to a lesser extent), exhibiting positive and slightly negative values along axis 2, respectively.

Mean values for the top 5 cm of the soil can be calculated for each class and each profile and averaged over the five profiles taken at each plot, giving mean values for the five plots (Table 2). The correspondence analysis displayed some trends common to gross classes (leaf material, charred material, charcoal, black and dark humus). These were compared between the five plots (Table 3). It appears that leaf material was at least 4 times more abundant in woody sites (MF and EF) than in abattis (MA, EA, EFB), each of these groups being fairly homogeneous. In particular, the topsoil of older abattis (EA, MA) had the same leaf litter content than the recently burnt forest (EFB). Slash-and-burn agriculture increased the content of the topsoil in black and dark hemorganic humus by a factor of 8 at Maripasoula. This class was 11 times more abundant at Elahe (EF) than at Maripasoula (MF) woody sites and remained 6 times more abundant at Elahe (EA) than in the abattis at Maripasoula (MA). Differences in the charcoal content of the topsoil were not so pronounced. Charcoal was less abundant in EF than in all other plots, the highest value being observed in MA (almost 2 times the value observed in EA).

Charred material was most abundant in the recently burnt forest (EFB), where it was 2.6 times more abundant than in the nearby untouched forest (EF). The least amount of charred material was in the Maripasoula woody fallow (MF).

The mean vertical distribution of these bulk classes is apparent in Figure 3. The leaf material, which was more abundant in woody sites (EF and MF) than in abattis, disappeared rapidly in the top 2 cm of the soil in all plots. The charred material, which was abundant at the soil surface in the two abattis at Elahe (EA, EFB), decreased with depth in these two plots. However, it disappeared much less rapidly than did the leaf material. In the other three plots, it increased progressively down to 1.5 to 2 cm, then decreased progressively. The vertical distribution of charcoal was similar to that of charred material in the recently burnt forest and in the cultivated abattis at Elahe, decreasing abruptly from 1.5 to 2.5 cm, then more progressively below this depth. In the Elahe forest (EF), there was a progressive increase in charcoal content from the surface to 2-3 cm, but the content remained very small. At Maripasoula (MA and MF), more charcoal had accumulated below 2.5 cm. In the abattis (MA), it reached on average 10% of the total solid matter at 2.5 cm, then decreased progressively.

In the Maripasoula forest (MF) the charcoal content increased progressively, from 0 at the surface to a maximum of 6% from 3 to 5 cm. Over this depth range, there was no difference in charcoal content between the 8-year-old woody fallow and the abattis. Black and dark humus was much more abundant in Elahe abattis (EA, EFB) than in all other plots. It was already present as a large proportion of total soil volume near the soil surface, reaching 9% and 6% of the total solid matter in EA and EFB, respectively. The proportion remained constant to 2 cm, then decreased progressively with increasing depth.

Discussion

Pronounced differences according to agricultural use of the rain forest were revealed by the biological material (leaf litter, charcoal and charred material, dark humus) which had accumulated in the topsoil. Biological activity in the soil was strongly affected by clearing of the forest. We registered a change from the activity of arthropods other than ants, as exemplified by the content in the topsoil of the faeces of caterpillars, millipedes, springtails and mites, to a community dominated by annelids and ants. The latter is characterized by (i) the incorporation of organic matter to mineral matter within earthworm and enchytraeid faeces and ant pellets, (ii) the building of an hemorganic humus profile (Figure 1). This pattern accords with the more even distribution of carbon in the topsoil after deforestation and cultivation observed in Amazonian soils by [START_REF] Nascimento | Evolution patterns of the soil organic matter in some agricultural systems in the Brazilian Cerrado region[END_REF]. Our results point on the role of burrowing animals, namely earthworms, enchytraeids and ants, in this process. We did not record any faecal deposition by termites, although many termites are present in the soil. The incorporation of organic matter in mineral matter by tropical earthworms helps to stabilize organic matter, and associated water and nutrients [START_REF] Lavelle | Large-scale effects of earthworms on soil organic matter and nutrient dynamics[END_REF]. Mutualistic interactions between earthworms and soil microbes (in particular the priming effect by intestinal mucus) were reasons to classify these animals as ecosystem engineers [START_REF] Lavelle | Soil function in a changing world: the role of invertebrate ecosystem engineers[END_REF].

Earthworms, in association with additions of organic matter, have been used for the bioremediation of tropical agricultural soils, increasing plant productivity and creating a stable soil structure [START_REF] Pashanasi | Effect of the endogeic earthworm Pontoscolex corethrurus on soil chemical characteristics and plant growth in a low-input tropical agroecosystem[END_REF][START_REF] Brown | Effects of earthworms on plant production in the tropics[END_REF][START_REF] Hallaire | Soil structure changes induced by the tropical earthworm Pontoscolex corethrurus and organic inputs in a Peruvian ultisol[END_REF]. Ants, and moreover enchytraeids, are not widely recognized as promoting greater soil aggregation and stability, at least in the tropical world [START_REF] Delabie | Soil and litter cryptic ant assemblages of Bahian cocoa plantations[END_REF][START_REF] Römbke | Applied research on Enchytraeidae in Central Amazonia: project approach, methodology and first results[END_REF]. However, we found that earthworm faeces tunnelled by enchytraeids contributed noticeably to the total solid matter in the first 5 cm of the topsoil (Table 2).

Charcoal and charred material contributed substantially to the soil (Figure 3). This material accumulates in the topsoil, unlike leaf material (Figure 3). It is also incorporated into a variety of black and dark hemorganic faeces of earthworms, and into ant pellets (Table 2). In a previous paper [START_REF] Topoliantz | Burrowing activity of the geophagous earthworm Pontoscolex corethrurus (Oligochaeta: Glossoscolecidae) in the presence of charcoal[END_REF], we showed experimentally that the peregrine earthworm Pontoscolex corethrurus could ingest small particles of charcoal and mix them with the mineral soil, then deposit dark excrement at varying depths. This activity is at the origin of a variety of dark colours of soil aggregates, which contain charcoal powdered by the earthworm gizzard then mixed in a liquid manner in the intestine [START_REF] Barois | Transformation of the soil structure through Pontoscolex corethrurus (Oligochaeta) intestinal tract[END_REF]. We also observed that ants produced dark, even black, material (Table 2), but this material could come from ingested earthworm faeces of the same colour.

There was more black and dark humus and charred material and less charcoal in the 3-year-old abattis at Elahe than in the 1-year-old abattis at Maripasoula. These differences can be due (i) to a larger production of charcoal at Maripasoula where repeated burning (short fallow) leads to an accumulation of charcoal in the course of the time ( [START_REF] Van Der Wal | Chinantec Shifting Cultivation: Interactive Landuse. A Case Study in the Chinantla, Mexico, on the Secondary Vegetation, Soils and Crop Performance under Indigenous Shifting Cultivation[END_REF], and (ii) to a better incorporation of charcoal and charred material in the mineral soil at Elahe, as shown by our micromorphological data. Indeed, charred material and charcoal decreased from EFB (fire 6 months before) to EA (fire 3 years before) while black and dark humus increased during the same period. This suggests that both charcoal and charred material are sources of black and dark humus (Figure 3). At Maripasoula, the content of the black and dark humus in the topsoil was always small and charcoal accumulated, revealing the slow turnover of this type of carbon.

The origin of the fertile Terra Preta, or Amazonian Dark Earths, has been attributed to the incorporation of charcoal to the mineral soil, as a source of stable humus, once charcoal carbon becomes oxidized and chemically reorganized [START_REF] Glaser | Black carbon in soils: the use of benzenecarboxylic acids as specific markers[END_REF]. Past bioturbation, in vanished civilizations of the Amazon Basin, has been postulated to explain the present-day enrichment of the soil in black carbon of charcoal origin at 30-40 cm depth [START_REF] Glaser | Black carbon in density fractions of anthropogenic soils of the Brazilian Amazon region[END_REF]Glaser & Woods, 2004). We can now explain this phenomenon by the charcoal-feeding activity of earthworms such as P. corethrurus, which are abundant in Amazonian pastures and forest clearings [START_REF] Römbke | About earthworm communities in a rain forest and an adjacent pasture in Peru[END_REF]. The fact that black and dark hemorganic material was almost absent from the soil in the old forest, then increased in the recently burnt forest, then still further increased after 3 years of cultivation at Elahe (Figure 3), probably indicates that the process of incorporation is rapid, provided earthworms are sufficiently active, as observed in Elahe. In previous experiments we also showed that dark hemorganic humus (artificially made by mixing charcoal with the mineral soil) was even ingested to a greater extent that pure charcoal, pointing on a positive feed-back (self-reinforcing) effect (Topoliantz & Ponge, 2005).

Our study shows that traditional shifting cultivation, as currently practised by Wayana Amerindians, causes charcoal and charred material to become well incorporated into the soil and to change gradually into stable humus. This probably occurred in the past to a much larger extent in the Amazonian basin (Glaser & Woods, 2004). Prior to burning the forest, Wayana cultivators select darker soils (personal communication with local agriculturists). The dark colour of the topsoil is probably caused by past fires that enriched the soil with stable carbon, through earthworm activity. This preference for soils containing incorporated charcoal reinforces the positive feed-back loop already mentioned.

At Maripasoula, the accumulation of charcoal in depth and the quasi absence of black and dark humus in the 8-years-old forest (MF) suggest that more than 8 years are probably needed for complete incorporation of charcoal in dark humus. Frequent fires and short cultivation periods must not allow for a complete incorporation of charcoal by soil fauna, limiting to the short-term the fertilizing effect of burning the forest. indebted to the Laboratory of Soil Analysis (INRA, Arras) for soil physicochemical analyses, and to various local people for commodities, free access to sampling sites, and fruitful discussions about agricultural practices and crops. C lay / % 2 2 . 4 (2 0 . 4 -2 5 . 3 ) 1 6 . 9 (1 3 . 7 -2 0 ) 2 0 (6 . 6 -3 1 . 3 ) 1 7 . 5 (1 5 . 5 -1 9 . 1 ) 1 6 . 8 (1 4 . 9 -1 8 . 3 ) S ilt / % 1 2 . 4 (9 . 5 -1 7 . 4 ) 7 . 9 (7 -9 . 2 ) 2 0 . 9 (1 3 . 2 -3 2 . 8 ) 8 . 7 (7 . 1 -9 . 6 ) 9 . 2 (7 . 6 -1 0 . 3 ) S an d / % 6 5 . 1 (5 8 . 8 -6 9 . 7 ) 7 5 . 2 (7 2 . 3 -7 8 . 7 ) 5 9 . 1 (3 5 . 9 -8 0 . 2 ) 7 3 . 8 (7 1 . 4 -7 6 . 7 ) 7 4 (7 2 . 1 -7 6 . 4 ) pH 5 . 1 3 (5 . 0 2 -5 . 3 7 ) 5 . 0 1 (4 . 6 -5 . 6 6 ) 4 . 7 4 (4 . 3 7 -4 . 9 1 ) 4 . 8 3 (4 . 5 8 -5 ) 4 . 3 0 (3 . 7 9 -5 . 5 9 )

T o t al C / g . k g -1 2 4 . 6 (1 6 . 9 -3 9 . 8 ) 2 2 . 5 (1 8 . 6 -2 8 . 4 ) 2 5 . 8 (1 6 . 3 -3 6 . 3 ) 1 9 . 1 (1 4 . 7 -2 5 . 3 ) 1 8 . 6 (1 7 . 9 -1 9 . 7 )

T o t al N / g . k g -1 1 . 6 5 (1 . 2 1 -2 . 4 9 ) 1 . 4 8 (1 . 2 8 -1 . 8 8 ) 1 . 7 8 (1 . 2 9 -2 . 4 ) 1 . 3 9 (1 . 1 2 -1 . 7 ) 1 . 3 5 (1 . 3 3 -1 . 3 8 )
C / N 1 4 . 7 (1 3 . 2 -1 6 ) 1 5 . 2 (1 4 . 5 -1 5 . 7 ) 1 4 . 6 (1 2 . 6 -1 6 . 1 ) 1 3 . 7 (1 2 . 8 -1 4 . 9 ) 1 3 . 7 (1 3 .

  , thus giving us confidence in the significance of these two axes. Axis 1 separated classes typical of woody sites (EF and MF, with positive values) from those typical of clearings (EA, EFB and MA, with negative values). Axis 2 separated EA (with positive values) from MA (with negative values), EFB being intermediate. The woody sites (EF and MF) were characterized by moss, leaf and wood litter at varying stages of decomposition, as well as a variety of holorganic (purely organic) faeces of arthropod origin (caterpillars, millipedes, springtails, mites, miscellaneous). Abattis (EA, EFB, MA) were characterized by annelid (earthworm, enchytraeid) faeces and ant pellets with varying carbon contents, charcoal (free or incorporated to faeces), roots and mineral particles. Within this group, the Elahe abattis (EA) was distinguished by the part played by dark (black, dark brown, dark grey) hemorganic material, and charred roots. Changes in the vertical distribution of topsoil components were shown by the projection of depth indicators in the plane of Axes 1 and 2 of CA (Figure 2).

Figure headings Figure 1 .

 headings1 Figure headings

Figure 2 . 1 Figure 3 .M

 213 Figure 2. Correspondence analysis (CA). Projection of depth indicators (for each plot and each depth) in the plane of the first two factorial axes. Successive sampling depths are shown at cm intervals. The top cm has been labelled, the lower samples follow the lines. Abbreviations as for Figure 1

  Figure 1 536 537

Table 2 .

 2 M ain p h y sic o -c h em ic al f eat u res o f t h e f iv e st u d y sit es. D at a are m ean s o f f iv e t o p so ils f o llo w ed b y m in im u m an d m ax im u m v alu es. C h arc o al c arb o n w as in c lu d ed in t o t al C List of classes identified under the dissecting microscope w ith code number and mean percent volume of the topsoil (0-5 cm) matrix in the five studied sites (5 profiles each). MA = Maripasoula Abattis, MF = Maripasoula Fallow , EA = Elahe Abattis, EF = Elahe Forest, EFB = Elahe Forest Burnt

	3 -1 4 . 3 )

Acknowledgements

Our study was financially supported by the Ministère de l'Écologie et du Développpement Durable (SOFT programme) and the Mission pour la Création du Parc de la Guyane, which are greatly acknowledged. We are greatly