Keywords: Problèmes hyperboliques du second ordre, méthodes explicites, méthodes réversibles en temps, conservation d'énergie, méthodes de Galerkine discontinu, méthodes d

Schéma à pas de temps local explicite et conservatif pour les équations d'ondes

Résumé : Les maillages raffinés localement imposent des fortes contraintes de stabilité aux méthodes explicites pour la résolution de problèmes hyperboliques. Pour surmonter cette difficulté, nous développons des méthodes de résolution permettant d'utiliser des petits pas de temps uniquement sur les éléments fins du maillage. Associée à une discrétisation en espace par éléments finis permettant d'obtenir une matrice de masse diagonale, cette méthode conduit à un schéma numérique explicite, parallélisable, et qui conserve une énergie discrète. À partir du schéma "saute-mouton" classique, nous construisons des méthodes à pas de temps local d'ordre arbitrairement élevé. Nous mettons en évidence l'utilité et l'efficacité de ces méthodes et nous validons nos résultats numériques au travers de résultats numériques.

Introduction

The efficient and accurate numerical solution of the wave equation is of fundamental importance for the simulation of time dependent acoustic, electromagnetic, or elastic wave phenomena. Finite difference methods are commonly used for the simulation of time dependent waves because of their simplicity and their efficiency on structured cartesian meshes [START_REF] Gustafsson | Time compact high order difference methods for wave propagation[END_REF][START_REF] Gustafsson | Time compact high order difference methods for wave propagation, 2[END_REF][START_REF] Kreiss | Difference approximations for the second order wave equation[END_REF]. However, in the presence of complex geometry or small geometric features that require locally refined meshes, their usefulness is somewhat limited. In contrast, finite element methods (FEMs) easily handle locally refined unstructured meshes; moreover, their extension to high order is straightforward, even in the presence of curved boundaries or material interfaces.

The finite element Galerkin discretization of second-order hyperbolic problems typically leads to a second-order system of ordinary differential equations. Even if explicit time-stepping is employed, the mass matrix arising from the spatial discretization by standard continuous finite elements must be inverted at each time-step: a major drawback in terms of efficiency. To overcome that problem, various "mass lumping" techniques have been proposed, which effectively replace the mass matrix by a diagonal approximation. While straightforward for piecewise linear elements [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], mass lumping techniques require special quadrature rules at higher order to preserve the accuracy and guarantee numerical stability [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF].

Alternatively, discontinuous Galerkin (DG) methods offer even greater flexibility for local mesh refinement by accomodating non-conforming grids and hanging nodes. Based on discontinuous finite element spaces, DG-FEMs weakly enforce continuity by adding suitable bilinear forms, so-called numerical fluxes, to standard variational formulations -see [START_REF] Cockburn | Discontinuous Galerkin methods for convection-dominated problems[END_REF][START_REF] Cockburn | The development of discontinuous Galerkin methods[END_REF][START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF] for further details and recent reviews. Because individual elements decouple, the mass matrix arising from a spatial DG discretization is block-diagonal, with block size equal to the number of degrees of freedom per element; it can therefore be inverted at very low computational cost. In fact, for a judicious choice of (locally orthogonal) shape functions, the mass matrix is diagonal. Thus, when combined with explicit time integration, the resulting time marching scheme will be truly explicit.

Recently, Grote, Schneebeli and Schötzau [START_REF] Grote | Discontinuous Galerkin finite element method for the wave equation[END_REF] proposed the symmetric interior penalty (IP) DG method for the second-order wave equation; in particular, they derived optimal a priori error bounds in the energy norm and the L 2 -norm for the semi-discrete formulation. A symmetric DG formulation of the wave equation in its second-order form offers the following advantage, which also pertains to the classical continuous Galerkin formulation. Since the stiffness matrix is positive semidefinite, the semi-discrete formulation inherently conserves (a discrete version of) the energy for all time. Moreover, when combined with a symmetric time marching scheme, such as the standard leap-frog (or Störmer-Verlet) method, the resulting fully discrete formulation will also conserve a discrete energy. Thus, both formulations will be free of any (unnecessary) damping. The dispersive properties of the symmetric IP-DG method were analyzed by Ainsworth, Monk and Muniz [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation[END_REF]. In [START_REF] Grote | Interior penalty discontinuous Galerkin method for Maxwell's equations: energy norm error estimates[END_REF][START_REF] Grote | Interior penalty discontinuous Galerkin method for Maxwell's equations: optimal L 2 -norm error estimates[END_REF] the symmetric interior DG method was extended to Maxwell's Equations in second-order form.

Adaptivity and mesh refinement are certainly key for the efficient numerical solution of partial differential equations. However, locally refined meshes impose severe stability constraints on explicit time-stepping schemes, where the maximal time-step allowed by the CFL condition is dictated by the smallest elements in the mesh. When mesh refinement is restricted to a small region, the use of implicit methods, or a very small time-step in the entire computational domain, are very high a price to pay. To overcome that stability restriction, various local time-stepping schemes were proposed, which allow smaller time-steps precisely where the smallest elements in the mesh are located. In [START_REF] Collino | A conservative space-time mesh refinement method for the 1-D wave equation[END_REF] Collino, Fouquet and Joly proposed a local time-stepping method for the first order wave equation; it was analyzed in [START_REF] Collino | A conservative space-time mesh refinement method for the 1-D wave equation[END_REF][START_REF] Joly | An error analysis of conservative space-time mesh refinement methods for the one-dimensional wave equation[END_REF] and extended to elastodynamics [START_REF] Bécache | Space-time mesh refinement for elastodynamics[END_REF] and Maxwell's equations [START_REF] Collino | Conservative space-time mesh refinement methods for the FDTD solution of Maxwell's equations[END_REF]. Their approach, which is based on the introdution of a Lagrange multiplier, conserves a discrete energy. However, it requires the solution of a linear system on the interface between the coarse and the fine mesh. By combining a symplectic integrator with a DG discretization of Maxwell's equations in first-order form, Piperno [START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problems[END_REF] proposed an explicit local time-stepping scheme, which also conserves a discrete energy. All these methods are second-order accurate in time. Alternatively, domain decomposition methods permit the use of different numerical methods or time steps in separate subdomains [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF][START_REF] Halpern | Local space-time refinement for the one-dimensional wave equation[END_REF].

The outline of our paper is as follows. In Section 2 we recall both the standard continuous and the symmetric IP-DG finite element discretizations of the scalar second-order wave equation, which typically serves as a model problem for general second-order hyperbolic problems. Starting from the well-known "leap-frog" scheme, we then derive a second-order a local timestepping scheme in Section 3. With a symmetric finite element discretization in space with a (block-)diagonal mass matrix, the resulting fully discrete scheme is not only explicit and thus inherently parallel, but it also conserves (a discrete version of) the energy. We also show via numerical experiments how a small overlap between the fine and the coarse region achieves an optimal CFL condition. Next, in Section 4, we extend the second-order local time-stepping schemes first to fourth order, and then to arbitrarily high order of accuracy. In Section 5 we present numerical experiments in one and two space dimensions which validate the theory and illustrate the usefulness of these local time-stepping schemes.

Finite element discretizations of the wave equation

We consider the scalar wave equation

u tt -∇ • c 2 ∇u = f in (0, T ) × Ω, (1)
u = 0 on (0, T ) × ∂Ω, (2)
u| t=0 = u 0 in Ω, (3)
u t | t=0 = v 0 in Ω, ( 4
)
where Ω is a bounded domain in R 2 or R 3 . Here f ∈ L 2 (0, T ; L 2 (Ω)) is a (known) source term, while u 0 ∈ H 1 0 (Ω) and v 0 ∈ L 2 (Ω) are prescribed initial conditions. We consider homogeneous Dirichlet boundary conditions, for simplicity, and assume that the speed of propagation, c(x), is piecewise smooth and strictly positive. In the absence of forcing, the (continuous) energy,

E[u](t) = 1 2 u t (t, .) 2 + ∇u(t, .) 2 , t ≥ 0, (5)
is conserved for all time.

We shall now discretize in space ( 1)-( 4) by using either standard continuous (H 1 -conforming) finite elements (with mass lumping) or a symmetric IP discontinuous Galerkin discretization from [START_REF] Grote | Discontinuous Galerkin finite element method for the wave equation[END_REF], while leaving time continuous. Thus, we consider shape-regular meshes T h that partition the domain Ω into disjoint elements {K}, such that Ω = ∪ K∈T h K. The elements are triangles or quadrilaterals in two space dimensions, and tetrahedra or hexaedra in three dimensions, respectively. The diameter of element K is denoted by h K and the mesh size, h, is given by h = max K∈T h h K .

Standard continuous Galerkin formulation

The standard continuous (H 1 -conforming) Galerkin formulation of the wave equation ( 1)-( 4) starts from its weak formulation [START_REF] Lions | Non-homogeneous Boundary Value Problems and Applications[END_REF]

: find u : [0, T ] × H 1 0 (Ω) → R such that (u tt , v) + (c∇u, c∇v) = (f, v) ∀v ∈ H 1 0 (Ω), t ∈ (0, T ), ( 6
)
u| t=0 = u 0 , (7)
u t | t=0 = v 0 . (8)
Here (., .) denotes the standard L 2 inner product on Ω. For a given partition T h of Ω and an approximation order ≥ 1, we wish to approximate the solution u(t, •) of ( 6)-( 8) in the finite element space

V h := {v ∈ H 1 0 (Ω) : v| K • F K ∈ S ( K) ∀K ∈ T h }, (9)
where S ( K) is the space P ( K) (for triangles or tetrahedra) or Q ( K) (for quadrilaterals or hexahedra) and F K : K → K is one-to-one and maps the boundary of the reference element, K, to the boundary of K. Thus, we consider the following semi-discrete Galerkin approximation of ( 6)-( 8): find

u h : [0, T ] × V h → R such that (u h tt , v) + (c∇u h , c∇v) = (f, v) ∀v ∈ V h , t ∈ (0, T ), (10)
u h | t=0 = Π h u 0 , (11)
u h t | t=0 = Π h v 0 , (12)
where Π h denotes the L 2 -projection onto V h -see [START_REF] Baker | Error estimates for finite element methods for second order hyperbolic equations[END_REF] for further details. Let (φ i ) i=1..N denote the standard nodal basis of V h defined by φ i (x j ) = δ ij , where (x i ) i=1..N are the nodes of the mesh. Next, we denote by y(t), F (t), y 0 , and ẏ0 the N -vectors

y i (t) = (u h (t, .), φ i ), F i (t) = (f (t, .), φ i ), y 0,i = (Π h u 0 , φ i ), ẏ0,i = (Π h v 0 , φ i ).
Then (10)-( 12) is equivalent to the second-order system of ordinary differential equations

M d 2 y dt 2 + Ky = F (13)
y(0) = y 0 , dy dt (0) = ẏ0 , (14)
where the N × N mass and stiffness matrices, M and K, are defined by

M ij = (φ i , φ j ), K ij = (c∇φ i , c∇φ j )),
respectively. The matrix M is symmetric positive definite whereas the matrix K is symmetric and, in general, positive semidefinite only.

Since the matrix M is not diagonal, it must be inverted at every time-step of any explicit time integration scheme. To overcome this difficulty, various so-called mass lumping techniques have been developed, which essentially replace M with a diagonal approximation by computing the integrals over each element K with judicious quadrature rules. For P 1 -elements in two space dimensions, for instance, one uses the approximation

K f |K| 3 3 k=1 f (x k K ),
where x k K are the vertices of triangle K. From the definition of the basis functions φ i , it then immediately follows that

K φ i φ j |K| 3 δ ij .
Hence, M is now diagonal while the spatial discretization remains second-order accurate [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF].

For P -elements up to order = 3, mass lumping techniques are also available but more complicated [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF], whereas the case of Q (or spectral)-elements for quadrilateral meshes and arbitrary , is well-understood [START_REF] Gary | Higher-order numerical methods for transient wave equations[END_REF][START_REF] Hennart | Topics in finite element discretization of parabolic evolution problems[END_REF][START_REF] Maday | Spectral element methods for the incompressible Navier-Stokes equations[END_REF][START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF]. Here, the degrees of freedom associated with the basis functions (φ i ) coincide with the Gauss-Lobatto quadrature points on each elements. As the integrals are also computed elementwise through Gauss-Lobatto quadrature, the resulting mass matrix is diagonal while the spatial accuracy is not affected [START_REF] Baker | Error estimates for finite element methods for second order hyperbolic equations[END_REF][START_REF] Baker | The effect of quadrature errors on finite element approximations for second order hyperbolic equations[END_REF].

Discontinuous Galerkin formulation

Here we briefly recall the symmetric interior penalty DG fomulation from [START_REF] Grote | Discontinuous Galerkin finite element method for the wave equation[END_REF]. For simplicity, we assume in this section that the elements are triangles or parallelograms in two space dimensions, and tetrahedra or parallelepipeds in three dimensions, respectively. Generally, we allow for irregular meshes with hanging nodes. An interior face of T h is the (nonempty) interior of ∂K + ∩ ∂K -, where K + and K -are two adjacent elements of T h . Similarly, a boundary face of T h is the (nonempty) interior of ∂K ∩ ∂Ω, which consists of entire faces of ∂K. We denote by F I h the set of all interior faces of T h , by F B h the set of all boundary faces, and let

F h = F I h ∪ F B h .
Here we generically refer to any element of F h as a "face", both in two and in three dimensions.

For any piecewise smooth function v we now introduce the following trace operators. Let F ∈ F I h be an interior face shared by two neighboring elements K + and K -and let x ∈ F ; we write n ± to denote the outward unit normal vectors on the boundaries ∂K ± . Denoting by v ± the trace of v taken from within K ± , we define the jump and average of v at x ∈ F by

[[v]] := v + n + + v -n -, { {v} } := (v + + v -)/2,
respectively. On every boundary face F ∈ F B h , we set [[v]] := vn and { {v} } := v. Here, n is the outward unit normal on ∂Ω.

For a piecewise smooth vector-valued function q, we analogously define the average across interior faces by { {q} } := (q + + q -)/2, and on boundary faces we set { {q} } := q. The jump of a vector-valued function will not be used. For a vector-valued function q with continuous normal components across a face f , the trace identity

v + (n + • q + ) + v -(n -• q -) = [[v]] • { {q} } on f,
immediately follows from the above definitions.

For a given partition T h of Ω and an approximation order ≥ 1, we wish to approximate the solution u(t, •) of ( 1)-(4) in the finite element space

V h := {v ∈ L 2 (Ω) : v| K ∈ S (K) ∀K ∈ T h }, (15)
where S (K) is the space P (K) of polynomials of total degree at most on K, if K is a triangle or a tetrahedra, or the space Q (K) of polynomials of degree at most in each variable on K, if K is a parallelogram or a parallelepiped. Thus, we consider the following (semi-discrete) discontinuous Galerkin approximation of (1)-( 4): find

u h : J × V h → R such that (u h tt , v) + a h (u h , v) = (f, v) ∀v ∈ V h , t ∈ J, (16)
u h | t=0 = Π h u 0 , (17)
u h t | t=0 = Π h v 0 . (18)
Here, Π h denotes the L 2 -projection onto V h and the discrete bilinear form a h on V h × V h is given by

a h (u, v) := K∈T h K c 2 ∇u • ∇v dx - F ∈F h F [[u]] • { {c 2 ∇v} } dA - F ∈F h F [[v]] • { {c 2 ∇u} } dA + F ∈F h F a [[u]] • [[v]] dA. ( 19
)
The last three terms in [START_REF] Fauqueux | Eléments finis mixtes spectraux et couches absorbantes parfaitement adaptées pour la propagation d'ondes élastiques en régime transitoire[END_REF] correspond to jump and flux terms at element boundaries; they vanish when u, v ∈ H 1 0 (Ω) ∩ H 1+σ (Ω), for σ > 1 2 . Hence the above semi-discrete discontinuous Galerkin formulation ( 16) is consistent with the original continuous problem [START_REF] Bécache | Space-time mesh refinement for elastodynamics[END_REF].

In [START_REF] Fauqueux | Eléments finis mixtes spectraux et couches absorbantes parfaitement adaptées pour la propagation d'ondes élastiques en régime transitoire[END_REF] the function a penalizes the jumps of u and v over the faces of T h . It is referred to as interior penalty stabilization function and is defined as follows. We first introduce the function h by

h| F = min{h K , h K }, F ∈ F I h , F = ∂K ∩ ∂K , h K , F ∈ F B h , F = ∂K ∩ ∂Ω. For x ∈ F , we further define c by c| F (x) = max{c| K (x), c| K (x)}, F ∈ F I h , F = ∂K ∩ ∂K , c| K (x), F ∈ F B h , F = ∂K ∩ ∂Ω.
Then, on each F ∈ F h , we set

a| F := α c 2 h -1 , (20)
where α is a positive parameter independent of the local mesh sizes and the coefficient c.

There exists a threshold value α min > 0 which only depends on the shape-regularity of the mesh and the approximation order , such that for α ≥ α min the discontinuous Galerkin bilinear form a h is coercive and, hence, the discretization stable [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF]. Throughout the rest of the paper we shall assume that α ≥ α min , so that the semi-discrete problem ( 16)-( 18) has a unique solution which converges with optimal order [START_REF] Grote | Discontinuous Galerkin finite element method for the wave equation[END_REF].

The semi-discrete IP-DG formulation ( 16)-( 18) is equivalent to the second-order system of ordinary differential equations

M d 2 y dt 2 + Ky = F (21)
y(0) = y 0 , dy dt (0) = ẏ0 . (22)
The N × N mass matrix M , with entries M ij = (φ i , φ j ), again is symmetric positive definite. Yet because individual elements decouple, M is also block-diagonal, with block size equal to the number of degrees of freedom per element. Thus, it can be inverted at very low computational cost. In fact, for a judicious choice of (locally orthogonal) shape functions, the mass matrix is diagonal and therefore the resulting time marching scheme truly explicit.

Remark 1. Because the matrices M and K are symmetic, either for the symmetric IP-DG or the standard continuous Galerkin discretization, the two semi-discrete formulations ( 13)-( 14) and ( 21)-( 22) conserve the (discrete) energy

E h (t) := 1 2 M y (t)y (t) + Ky(t), y(t) .
when F = 0. Here the angular brackets denote the standard Euclidean inner product on R N . If the underlying bilinear form a is strictly coercive, so are the two discrete bilinear forms a h , either for the discontinuous Galerkin discretization with α ≥ α min or for the continuous Galerkin discretization with (or without) mass lumping. For f = 0, conservation of energy then implies that the solution remains bounded for all time. However, when the underlying bilinear form a is not (strictly) coercive, as in the presence of Neumann or periodic boundary conditions, for instance, the elliptic partial differential operator will have a (simple) zero eigenvalue (with constant eigenfunction). Then, the wave equation admits linear growth in time, although the (continuous or discrete) energy remains constant. Nonetheless, if the numerical method is consistent, stability (with respect to the energy) will imply convergence on any finite time interval. Moreover, since all time integration schemes considered here are at least second-order accurate, the numerical time integration of the linearly growing zeroth eigenmode will actually be exact.

Local time-stepping

We consider the semi-discrete wave equation,

M d 2 y dt 2 + Ky = 0, ( 23
)
where M is an N × N symmetric positive definite (sparse) matrix and K is an N × N symmetric positive semidefinite (sparse) matrix. Moreover, we assume that M 1 2 can be explicitly computed and inverted at low cost, as in the case when M is (block-) diagonal, for instance. Next, we multiply ( 23) by M -1 2 to obtain

d 2 z dt 2 + M -1 2 KM -1 2 z = 0, (24)
with Let A denote the matrix M -1 2 KM -1 2 , which is also sparse and symmetric positive semidefinite. Then, we rewrite [START_REF] Grote | Interior penalty discontinuous Galerkin method for Maxwell's equations: optimal L 2 -norm error estimates[END_REF] as

d 2 z dt 2 + Az = 0 . ( 25
)
For any f ∈ C 1 we have

f (t + ∆t) -2f (t) + f (t -∆t) = -∆t 2 1 -1 (1 -|θ|)f (t + θ∆t) dθ. (26)
Hence the (exact) solution z(t) of ( 25) satisfies

z(t + ∆t) -2z(t) + z(t -∆t) = -∆t 2 1 -1 (1 -|θ|)Az(t + θ∆t) dθ. (27)

Second-order local time-stepping method

The integral on the right side of ( 27) represents a weighted average of Az(s) over the interval [t -∆t, t + ∆t], which needs to be approximated in any numerical algorithm. For instance, if we simply replace Az(t + θ∆t) by Az(t) in ( 27) and evaluate the remaining θ-dependent integral, we obtain the well-known second-order leap-frog scheme with time-step ∆t,

z n+1 -2z n + z n-1 = -∆t 2 Az n , z n z(t n ), (28)
which, however, would require ∆t to be comparable in size to the smallest elements in the mesh. Instead, we now split the vector z(t) in two parts:

z(t) = (I -P )z(t) + P z(t) = z [coarse] (t) + z [fine] (t) .
The projection matrix P is diagonal: its diagonal entries, equal to zero or one, identify the unknowns associated with the locally refined region, that is where smaller time-steps are needed. To circumvent the severe CFL restriction on ∆t in [START_REF] Halpern | Local space-time refinement for the one-dimensional wave equation[END_REF], we shall treat z

[fine] (t) differently from z [coarse] (t) in z(t + ∆t) -2z(t) + z(t -∆t) = -∆t 2 1 -1 (1 -|θ|) A z [coarse] (t + θ∆t) + z [fine] (t + θ∆t) dθ. (29)
Following [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Leimkuhler | Simulating Hamiltonian dynamics[END_REF], we approximate the integrand in [START_REF] Hennart | Topics in finite element discretization of parabolic evolution problems[END_REF] as

A z [coarse] (t + θ∆t) + z [fine] (t + θ∆t) Az [coarse] (t) + AP z(θ∆t),
where z(τ ) solves the differential equation

       d 2 z dτ 2 (τ ) = -A(I -P )z(t) -AP z(τ ), z(0) = z(t), dz dτ (0) = ν, (30)
and ν will be specified below. Thus,

z(t + ∆t) -2z(t) + z(t -∆t) -∆t 2 1 -1 (1 -|θ|) [A(I -P )z(t) + AP z(θ∆t)] dθ. (31)
Note that A and P do not commute. Since z solves [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], we deduce again from (26) that

z(∆t) -2z(0) + z(-∆t) = -∆t 2 1 -1 (1 -|θ|) [A(I -P )z(t) + AP z(θ∆t)] dθ. ( 32
)
From the comparison of ( 31) and [START_REF] Kreiss | Difference approximations for the second order wave equation[END_REF] we infer that

z(t + ∆t) -2z(t) + z(t -∆t) z(∆t) -2z(0) + z(-∆t). ( 33
)
The quantity z(∆t) -2z(0) + z(-∆t) does not depend on the value of ν, which we choose equal to zero. Therefore, z(∆t) = z(-∆t) by symmetry and (33) becomes

z(t + ∆t) -2z(t) + z(t -∆t) 2(z(∆t) -z(0)), (34)
or equivalently z(t + ∆t) + z(t -∆t) 2z(∆t).

Hence, we shall approximate the right side of (35) by solving (30) on [0, ∆t], and then use [START_REF] Lubich | [END_REF] to compute z(t + ∆t). Since the first term on the right side of (30) does not depend on τ , the (high-frequency, oscillatory) evolution of z(τ ) is solely determined by the second term, AP z(τ ), which involves only the unknowns associated with the nonzero entries in P . If those nonzero entries occupy only a small fraction of all unknowns, the additional effort from solving (30) will be small since A is sparse. Clearly, in doing so we must also ensure that the overall numerical scheme remains second-order accurate in time, as we shall show below. In summary, the local time-stepping algorithm for the solution of ( 23) computes z n+1 z(t + ∆t), given z n and z n-1 , as follows: Here Steps 1-3 correspond to the numerical solution of (30) with ν = 0 until τ = ∆t using the leap-frog scheme with the local time-step ∆τ = ∆t/p. In fact, any other second-order method, either explicit or implicit, could be used there instead. For P = 0, that is without any local time-stepping, we have

z1 = z(∆t) = I - ∆t 2 2 A z(t)
and hence we recover the standard leap-frog scheme [START_REF] Halpern | Local space-time refinement for the one-dimensional wave equation[END_REF]. If the fraction of nonzero entries in P is small, the overall cost will be dominated by the computation of w, which requires a single multiplication by A(I -P ) per time-step ∆t. All further matrix-vector multiplications by AP involve only those unknowns that are associated with the smaller, locally refined region. In addition, since A is sparse every update in Step 3 affects only those unknowns that lie inside the refined region, or immediately next to it.

To establish the accuracy and stability of the above local time-stepping scheme we shall now show how to rewrite it in "leap-frog manner". To do so, we first need the following technical result. Lemma 3. For m ≥ 2, zm/p defined by Algorithm 2 satisfies

zm/p = z n - m 2 2 ∆t p 2 Az n + m-1 j=1 ∆t p 2(j+1) α m j (AP ) j Az n , (37)
where the constants α m j are given by

                           α 2 1 = 1 2 , α 3 1 = 3, α 3 2 = - 1 2 , α m+1 1 = m 2 2 + 2α m 1 -α m-1 1 , α m+1 j = 2α m j -α m-1 j -α m j-1 , j = 2...m -2, α m+1 m-1 = 2α m m-1 -α m m-2 , α m+1 m = -α m m-1 . (38)
Proof. The proof is by induction on m. We first show that (37) holds for m = 2. Since

z1/p = z n - 1 2 ∆t p 2 Az n ,
we immediately find that

z2/p = z n -2 ∆t p 2 Az n + 1 2 ∆t p 4 AP Az n .
Hence, (37) holds with α 2 1 = 1/2. Next, let (37) hold for m. Then,

z(m+1)/p = 2   zn - m 2 2 ∆t p 2 Az n + m-1 j=1 ∆t p 2(j+1) α m j (AP ) j Az n    -   zn - (m -1) 2 2 ∆t p 2 Az n + m-2 j=1 ∆t p 2(j+1) α m-1 j (AP ) j Az n    - ∆t p 2   A(I -P )z n + AP   zn - m 2 2 ∆t p 2 Az n + m-1 j=1 ∆t p 2(j+1) α m j (AP ) j Az n       ,
which after some algebra simplifies to

z(m+1)/p = z n -m 2 - (m -1) 2 2 + 1 ∆t p 2 Az n + m 2 2 ∆t p 4 AP Az n + 2 m-1 j=1 ∆t p 2(j+1) α m j (AP ) j Az n - m-2 j=1 ∆t p 2(j+1) α m-1 j (AP ) j Az n - m-1 j=1 ∆t p 2(j+2) α m j (AP ) j+1 Az n .
Rearranging terms we then find

z(m+1)/p = z n - (m + 1) 2 2 ∆t p 2 Az n + ∆t p 4 m 2 2 + 2α m 1 -α m-1 1 AP Az n + m-2 j=2 ∆t p 2(j+1) 2α m j -α m-1 j -α m j-1 (AP ) j Az n + ∆t p 2m 2α m m-1 -α m m-2 (AP ) m-1 Az n - ∆t p 2(m+1) α m m-1 (AP ) m Az n ,
which yields [START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problems[END_REF] with α m j as in [START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF].

As a consequence, we can rewrite the above local time-stepping algorithm in "leap-frog manner". Proposition 4. The local time-stepping Algorithm 2 is equivalent to

z n+1 = 2z n -z n-1 -∆t 2 A p z n , (39)
where A p is defined by

A p = A - 2 p 2 p-1 j=1 ∆t p 2j α p j (AP ) j A (40)
and the constants α p j are given by [START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF]. This scheme is second-order accurate. Furthermore, the matrix A p is symmetric.

Proposition 4 is crucial for the accuracy and stability analysis below. However, the actual implementation of the local time-stepping scheme follows Algorithm 2; in particular, neither A p nor the constants α p j are ever used in practice.

Proof. Recall that z n+1 = -z n-1 + 2z 1 . We now use [START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problems[END_REF] with m = p to replace z1 . This yields

z n+1 = -z n-1 + 2z n -∆t 2 Az n + 2 ∆t p 2 p-1 j=1 ∆t p 2j α p j (AP ) j Az n , (41)
which corresponds to [START_REF] Shubin | A modified equation approach to constructing fourth-order methods for acoustic wave propagation[END_REF] with A p as in [START_REF] Villamizar | Time-dependent numerical method with boundary-conforming curvilinear coordinates applied to wave interactions with prototypical antennas[END_REF].

To prove consistency we now rewrite (41) as

z n+1 -2z n + z n-1 ∆t 2 - 2 p 2 p-1 j=1 ∆t p 2j α p j (AP ) j Az n = -Az n . Since d 2 z dt 2 (t) = z(t + ∆t) -2z(t) + z(t -∆t) ∆t 2 + O(∆t 2 ) and p-1 j=1 ∆t p 2j α p j (AP ) j Az(t) = O(∆t 2 ),
we conclude that the local time-stepping scheme is second-order accurate in time.

Finally, as the matrices A and P are symmetric, we have

((AP ) j A) = A(P A) j = (AP ) j A, j ≥ 1.
Therefore, the matrix A p is symmetric, too.

Energy conservation

The standard leap-frog scheme [START_REF] Halpern | Local space-time refinement for the one-dimensional wave equation[END_REF] conserves the discrete energy,

E n+ 1 2 = 1 2 I - ∆t 2 4 A z n+1 -z n ∆t , z n+1 -z n ∆t + A z n+1 + z n 2 , z n+1 + z n 2 . ( 42
)
Here

E n+ 1 2 E(t n+ 1 2
) (see Remark 3.2) and the angular brackets denote the standard Euclidean inner product on R N . Since A is symmetric, so is the quadratic form in (42). Moreover, for sufficiently small ∆t it is also positive semidefinite and hence yields a true energy. That restriction on ∆t corresponds precisely to the CFL stability condition of the leap-frog method and guarantees its numerical stability.

To prove the numerical stability of the local time-stepping algorithm, and thereby determine a necessary and sufficient condition for stability, we shall now exhibit a conserved discrete energy.

Proposition 5. The second-order local time-stepping scheme (Algorithm (2)) conserves the following energy:

E n+ 1 2 = 1 2 I - ∆t 2 4 A p z n+1 -z n ∆t , z n+1 -z n ∆t + A p z n+1 + z n 2 , z n+1 + z n 2 . ( 43
)
This energy corresponds to the energy conserved by the leap-frog scheme with A replaced by A p .

Proof. First, we take the inner product of (39) with z n+1 -z n-1 and divide the resulting expression by ∆t 2 to obtain

z n+1 -2z n + z n-1 ∆t , z n+1 -z n-1 ∆t + A p z n , z n+1 -z n-1 = 0.
By using

z n+1 -2z n + z n-1 = (z n+1 -z n ) -(z n -z n-1 )
and

z n+1 -z n-1 = (z n+1 -z n ) + (z n -z n-1 )
we obtain

z n+1 -z n ∆t , z n+1 -z n ∆t - z n -z n-1 ∆t , z n -z n-1 ∆t + A p z n , z n+1 -z n-1 = 0.
Since A p is symmetric, we find

A p z n , z n+1 -z n-1 = A p z n+1 , z n -A p z n , z n-1 ,
which proves that the quantity

E n+ 1 2 = z n+1 -z n ∆t , z n+1 -z n ∆t + A p z n+1 , z n is conserved. Next, we remark that A p z n+1 , z n = A p z n+1 + z n 2 , z n+1 + z n 2 - ∆t 2 4 A p z n+1 -z n ∆t , z n+1 -z n ∆t .
Thus, we have

E n+ 1 2 = I - ∆t 2 4 A p z n+1 -z n ∆t , z n+1 -z n ∆t + A p z n+1 + z n 2 , z n+1 + z n 2 .
The above local time-stepping algorithm conserves the energy E n+ For p = 1 we have A p = A, and thus we recover the well-known CFL condition of the standard leap-frog scheme:

∆t ≤ 2 √ λ max = ∆t LF .
For p > 1, the matrix A p explicitly depends on ∆t, and so do its eigenvalues. Moreover, as the eigenvectors of A and A p generally do not coincide, the analytic derivation of a CFL condition is not obvious. Instead, we shall perform a systematic numerical study of the eigenvalues of

∆t 2
4 A p in the following typical situation.

Stability and CFL condition: numerical study

We consider the one-dimensional wave equation with constant wave speed c = 1 on the interval Ω = [0 ; 6] with periodic boundary conditions. Next, we divide Ω into three equal parts. The left and right intervals, [0 ; 2] and [4 ; 6], are discretized with an equidistant mesh of size h coarse , whereas the interval [2 ; 4] is discretized with an equidistant mesh of size h fine = h coarse /p. Hence, the two outer intervals correspond to the coarse region and the inner interval, [2 ; 4], to the refined region -see Fig. 1.

For every time-step ∆t, we shall take p steps of size ∆τ = ∆t/p in the refined region. In the absence of local refinement, i.e. p = 1, the mesh is equidistant throughout Ω. Then, the (local) time-stepping algorithm corresponds to the standard leap-frog (LF) method and we denote by ∆t LF the largest time-step allowed. For p ≥ 2, we let ∆t p denote the maximal timestep of Algorithm 2. If ∆t p = ∆t LF , the local time-stepping algorithm imposes no further restriction on ∆t and we then shall call the CFL condition of the new scheme optimal.

First, we consider the IP-DG discretization with P 1 -elements and (small) penalization, α = 2. We choose h coarse = 0.2, which yields the maximal time-step ∆t LF = 0.55

h coarse = h coarse h fine 0 2 4 6 z [fine]
z [fine] with an overlap by one element Figure 1: IP-DG P 1 -elements. The computational mesh and the associated degrees of freedom, indicated by small arrows. 0.11 for p = 1. Now, we refine by a factor p = 2 those elements that lie inside the interval [START_REF] Anné | Construction and analysis of higher order finite difference schemes for the 1D wave equation[END_REF][START_REF] Baker | Error estimates for finite element methods for second order hyperbolic equations[END_REF], that is h fine = 0.1, and set to one all corresponding entries in P . Hence for every time-step ∆t, we shall take two steps of size ∆τ = ∆t/2 in the refined region.

To determine the range of values ∆t for which the local time-stepping scheme is stable, we display the eigenvalues of (∆t 2 /4)A p for varying ∆t/∆t LF -recall that A p also depends on ∆t. The numerical scheme is stable for any particular ∆t if all corresponding eigenvalues lie between zero and one; otherwise, it is unstable. Since the smallest eigenvalue actually never dips below zero, it does not affect the stability here. As shown in the left frame of Fig. 2, the largest time-step allowed is only about 60% of ∆t LF ; hence, the gain over a straightforward reduction of the (global) time-step by a factor two is rather modest.

To allow for larger time-steps, we now slightly enlarge the set of unknowns where a local time-step is used by also including those degrees of freedom that are associated with elements directly adjacent to the refined region. By setting the corresponding entries in P to one, we easily realize this overlap by one element in z [fine] . In the right frame of Fig. 2 we observe that all eigenvalues now lie essentially between zero and one. However, a hundred thousandfold magnification of that same figure, shown in the left frame of Fig. 3, reveals that some eigenvalues still barely transgress the strict stability limit at one. Further extension of the overlap by one additional element removes all unstable values below 0.9∆t LF , as shown in the right frame of Fig. 3, while four narrow bands of (barely) unstable values between 0.91 ≤ ∆t/∆t LF ≤ 0.98 remain. Here we shall not attempt to elucidate that peculiar and somewhat sensitive behavior, due to weak resonances caused by the underlying regular, one-dimensional grid. Instead, we now address the question whether the local time-stepping scheme is stable for the maximal time-step ∆t = ∆t LF . In Tab. 1 we list the corresponding maximal eigenvalue of (∆t 2 /4)A p for an overlap by one element with different mesh sizes h coarse = 0.5, 0.2, 0.1, 0.05, 0.025 and different p = 2, 3, 4, 10, 13; here, the maximal eigenvalue typically is greater than one and the local time-stepping scheme therefore unstable for the optimal The eigenvalues of (∆t 2 /4)A p : overlap by one element (left); overlap by two elements (right). The vertical scale is strongly magnified: 0.99999 < λ max < 1.00001. elements (IP-DG P 1 elements). time-step. Yet with an overlap by two elements, the maximal eigenvalues listed in Tab. 2 now always lie below one, independently of h coarse and p; hence for an overlap by two elements, Algorithm 2 is always stable for the optimal time-step. 0.9999 0.9999 0.9999 0.9999 0.9999 0.025 0.9999 0.9999 0.9999 0.9999 0.9999 Next, we repeat the above experiment with standard P 1 continuous, piecewise linear finite elements with mass lumping. Then, the CFL condition for the classical leap-frog scheme is ∆t LF = h coarse = 0.2. In contrast to the IP-DG FE method, a degree of freedom can now belong to both a fine and a coarse element; hence, any degree of freedom at the interface between the two sub-regions is automatically included in z [fine] , as shown in Fig. 4.

0 2 4 6 ¡ ¢ £ ¤ ¥ §¦ © § © § © § © § ! " # z [fine]
z [fine] with an overlap by one element Figure 4: The nodes in z [fine] for continuous P 1 elements

In Tab. 3 we list the maximal eigenvalue of ∆t 2 4 A p at ∆t = ∆t LF for an overlap by a single element for different mesh sizes h coarse = 0.5, 0.2, 0.1, 0.05, 0.025 and different p = 2, 3, 4, 10, 13. Here, we find that an overlap by only one element already leads to a stable scheme for the optimal time-step. In summary, a slight extension (overlap) of the region where local time-steps are used into that part of the mesh immediately adjacent to the refined region always improves the stability of the time-stepping scheme. An overlap by one element for the P 1 continuous FE discretization (with mass lumping), or by two elements for the IP-DG discretization, always permits the use of the maximal (optimal) time-step ∆t LF , dictated by the leap-frog method in the coarse region.

High order local time-stepping

We shall now show how to extend the second-order local time-stepping method from Section 3 to arbitrarily high accuracy. First, we develop in detail a fourth-order local time-stepping scheme and again exhibit a conserved discrete energy. Then, we consider the general case of arbitrary (even) order. By including a small overlap into the surrounding coarser region we show via numerical experiments that both the fourth-and sixth-order local time-stepping schemes are stable for the optimal time-step dictated by the coarse mesh size.

Fourth-order local time-stepping

A fourth-order extension of the leap-frog scheme for ( 25) is provided by the Modified Equation (ME) approach [START_REF] Carpentier | On the derivation of the modified equation for the analysis of linear numerical methods[END_REF][START_REF] Dablain | High order differencing for the scalar wave equation[END_REF][START_REF] Anné | Construction and analysis of higher order finite difference schemes for the 1D wave equation[END_REF][START_REF] Shubin | A modified equation approach to constructing fourth-order methods for acoustic wave propagation[END_REF]. Thus, we replace Az(t + θ∆t) in ( 27) by the leading terms in its Taylor expansion,

Az(t + θ∆t) A z(t) + θ∆tz (t) + θ 2 ∆t 2 2 z (t) .
By integrating in [START_REF] Hairer | Geometric Numerical Integration[END_REF] the resulting polynomial in θ, the first derivative of z vanishes; since z (t) = -Az(t), we thus obtain the modified equation scheme:

z m+1 -2z m + z m-1 ∆t 2 = -Az m + ∆t 2 12 A 2 z m , z m z(t m ). ( 44
)
In contrast to standard fourth-order Runge-Kutta schemes, the ME approach requires only two multiplications with A per time-step. Its CFL condition follows from the CFL condition of the leap-frog scheme [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF],

∆t M E ≤ 2 √ 3 √ λ max = √ 3∆t LF .
As the maximal time-step allowed by the modified equation approach is about 70% times larger than that of the leap-frog scheme, the additional work needed for the improved accuracy is quite small. Clearly, high-order schemes of arbitrary (even) order can be obtained by using additional terms in the Taylor expansion of Az(t + θ∆t). Recently Gilbert and Joly derived new time-stepping methods that maximize the CFL condition for a given order of accuracy [START_REF] Gilbert | Higher order time stepping for second order hyperbolic problems and optimal cfl conditions[END_REF].

Following [START_REF] Lubich | [END_REF], we expand in (29) z [coarse] in Taylor series as

z [coarse] (t + θ∆t) = z [coarse] (t) + θ∆t dz [coarse] dt (t) + θ 2 ∆t 2 2 d 2 z [coarse] dt 2 (t) + O(∆t 3 ). Since d 2 z [coarse] dt 2 (t) = (I -P ) d 2 z dt 2 (t) = -(I -P )Az(t),
we find that

z [coarse] (t + θ∆t) (I -P )z(t) + θ∆t(I -P ) dz dt (t) - θ 2 ∆t 2 2 (I -P )Az(t). ( 45
)
Next, we use ( 45) in ( 29) to replace z [coarse] , which yields

z(t + ∆t) -2z(t) + z(t -∆t) -∆t 2 A(I -P )z(t) + ∆t 4 12 A(I -P )Az(t) -∆t 2 1 -1 (1 -|θ|)Az [fine] (t + θ∆t)dθ. ( 46
)
Hence if P = 0 we recover the standard modified equation scheme (44).

Similarly to Section 3.1, we shall now approximate the right-hand side of (46) by solving the following differential equation for z(τ ):

       d 2 z dτ 2 = -A(I -P )z(t) + τ 2 2 A(I -P )Az(t) -AP z(τ ) z(0) = z(t), dz dτ (0) = 0. ( 47
)
Again by symmetry, we have

z(t + ∆t) + z(t -∆t) 2z(∆t). (48)
Thus, we approximate the right side of (46) by solving (47) on [0, ∆t], and then use (48) to compute z(t + ∆t). The last term on the right of (46) explicitly depends on z(τ ), but it involves only the unknowns associated with the nonzero entries in P . If those nonzero entries occupy only a small fraction of all unknowns, the additional effort from solving (47) will be small. Clearly, in doing so we must ensure that the overall numerical scheme remains fourth-order accurate in time.

In summary, the fourth-order local time-stepping algorithm for the solution of ( 23) computes z n+1 z(t + ∆t), given z n and z n-1 , as follows:

Algorithm 6. 1. Set z0 = z n , w 1 = A(I -P )z n , w 2 = A(I -P )Az n , u 1 = AP z n and u 2 = AP Az n ; 2. Compute z1/p = z0 -1 2 ∆t p 2 (w 1 + u 1 ) + 1 24 ∆t p 4 (w 2 + u 2 ); 3. For m = 1...p -1, compute v 1 = -w 1 + 1 2 m∆t p 2 w 2 -AP zm/p ; v 2 = w 2 -AP v 1 ; z(m+1)/p = 2z m/p -z(m-1)/p + ∆t p 2 v 1 + 1 12 ∆t p 4 v 2 ; 4. Compute z n+1 = -z n-1 + 2z 1 .
Here Steps 1-3 correspond to the numerical solution of (47) until τ = ∆t with the modified equation approach using the local time-step ∆τ = ∆t/p. This algorithm requires two multiplications by A(I -P ) in Step 1 and 2p further multiplications by AP . No multiplication by A is needed to compute Az n in Step 1, since

Az n = A(I -P )z n + AP z n = w 1 + u 1 .
For P = 0, that is without any local time-stepping, this algorithm reduces to the modified equation scheme (44) above.

To establish the accuracy and stability of the above local time-stepping scheme we shall now show how to rewrite it in "leap-frog manner". To do so, we first need the following technical result. Lemma 7. For m ≥ 2, zm/p defined by Algorithm 6 satisfies

zm/p = z n - m 2 2 ∆t p 2 Az n + m 4 24 ∆t p 4 A 2 z n + 2(m-1) j=1 ∆t p 2(j+2) β m j (AP ) j AAz n (49)
where β m j are constant. Proof. The proof is by induction over m.

We first show that (49) holds for m = 2. Starting from the definitions of w 1 and w 2 in Step 1 of Algorithm 6, we find that

v 1 = -A(I -P )z n + 1 2 ∆t p 2 A(I -P )Az n -AP z1/p , v 2 = A(I -P )Az n -AP v 1 . Since z1/p = z n - 1 2 ∆t p 2 Az n + 1 24 ∆t p 4 A 2 z n ,
the two expressions for v 1 and v 2 simplify as

v 1 = -Az n + 1 2 ∆t p 2 A 2 z n - 1 24 ∆t p 4 AP A 2 z n , (50)
v 2 = A 2 z n - 1 2 ∆t p 2 AP A 2 z n + 1 24 ∆t p 4 AP AP A 2 z n . (51)
We now replace v 1 , v 2 in the expression for z2/p in Step 3. This yields

z2/p = z n -2 ∆t p 2 Az n + 2 3 ∆t p 4 A 2 z n (52) - 1 12 ∆t p 6 AP A 2 z n + 1 288 ∆t p 8 AP AP A 2 z n . (53)
Hence, β 2 1 = -1/12 and β 2 1 = 1/288, which completes the proof for m = 2. For higher m ≥ 3 (and p ≥ 3), the proof iso by induction on m; those straightforward but tedious calculations are omitted here.

As a consequence, we can rewrite the above local time-stepping algorithm in "leap-frog manner". Proposition 8. The local time-stepping Algorithm 6 is equivalent to

z n+1 = 2z n -z n-1 -∆t 2 A p z n , (54)
where A p is defined by

A p = A - ∆t 2 12 A 2 - 2 p 2 2(p-1) j=1 ∆t p 2(j+1) β p j (AP ) j A 2 . ( 55
)
This scheme is fourth-order accurate. Furthermore, the matrix AA p is symmetric.

Proof. Recall that z n+1 = -z n-1 + 2z 1 . We now use (49) with m = p to replace z1 . This yields

z n+1 = -z n-1 + 2z n -∆t 2 Az n + ∆t 4 12 A 2 z n + 2∆t 2 p 2 2(p-1) j=1 ∆t p 2(j+1) β p j (AP ) j AAz n ,
which corresponds to (54) with A p as in (55).

To prove consistency, we rewrite the scheme as

z n+1 -z n-1 + 2z n ∆t 2 - ∆t 2 12 A 2 z n - 2 p 2 2(p-1) j=1 ∆t p 2(j+1) β p j (AP ) j AAz n = -Az n . From z(t + ∆t) -2z(t) + z(t -∆t) ∆t 2 - ∆t 2 12 A 2 z(t) = z (t) + O(∆t 4 )
and

2(p-1) j=1 ∆t p 2(j+1) β p j (AP ) j A 2 z(t) = O(∆t 4 ),
we conclude that Algorithm 6 is indeed a fourth-order approximation of [START_REF] Gustafsson | Time compact high order difference methods for wave propagation[END_REF].

As the terms (AP ) j A 2 are obviously not symmetric, the matrix A p is not symmetric either. However, since (A(AP ) j A) = A(P A) j A = (AP ) j A, j ≥ 1, the matrix AA p indeed is symmetric.

From Tab. 4 we conclude that Algorithm 6 without overlap is unstable for the maximal timestep. However, as we extend by a single adjacent element the region where smaller time-steps are used, the fourth-order local time-stepping scheme becomes stable independently of p and h coarse .

h coarse AA p A -(∆t 2 /4)AA p 0.5 1.2 1.1 0.2 1.2 1.1 0.1 1.2 1.1 0.05 1.2 1.1 0.025 1.2 1.1
Next, we consider continuous piecewise P 3 finite elements with mass lumping. Since the CFL condition of the leap-frog scheme is ∆t LF = 0.232 h coarse [START_REF] Fauqueux | Eléments finis mixtes spectraux et couches absorbantes parfaitement adaptées pour la propagation d'ondes élastiques en régime transitoire[END_REF], the corresponding CFL condition of the modified equation scheme is ∆t M E = √ 3∆t LF = 0.401 h coarse . Again, we automatically include those degrees of freedom that belong to both a fine and a coarse cell into z [fine] . Remarkably, no overlap is needed here for the local time-stepping scheme to remain stable, regardless of and h coarse , as shown in Tab. 5 for p = 2.

Local time-stepping of arbitrary order

The modified equation approach used in Section 4.1 to extend the local time-stepping approach to fourth-order can be generalized to arbitrarily accuracy. Since the procedure parallels that used in Section 4.1, we shall omit details and state only key results.

Starting again from the Taylor expansion

z [coarse] (t + θ∆t) = z [coarse] (t) + 2s-2 i=1 θ i ∆t i i! d i z [coarse] dθ i (t) + O(∆t 2s-1 )),
we use the fact that

d 2i z [coarse] dt 2i (t) = (I -P ) d 2i z dt 2i (t) = (I -P )(-A) i z(t), to obtain z [coarse] (t + θ∆t) = z [coarse] (t) + s-1 i=1 θ 2i ∆t 2i (2i)! (I -P )(-A) i z(t) + s-1 i=1 θ 2i-1 ∆t 2i-1 (2i -1)! d 2i-1 z [coarse] dθ 2i-1 (t) + O(∆t 2s-1 )). (58)
Next, we use (58) in [START_REF] Hennart | Topics in finite element discretization of parabolic evolution problems[END_REF] to replace z [coarse] , which yields

z(t + ∆t) -2z(t) + z(t -∆t) -∆t 2 A(I -P )z(t) + 2 s-1 i=1 ∆t 2i+2 (2i + 2)! A(I -P )(-A) i z(t) -∆t 2 1 -1 (1 -|θ|)Az [fine] (t + θ∆t)dθ. (59)
Similarly to Section 3.1, we now approximate the right-hand side of (59) by solving

       d 2 z dτ 2 = -A(I -P )z(t) - s i=1 τ 2i 2i! A(I -P )(-A) i z(t) -AP z(τ ), z(0) = z(t), dz dt (0) = 0, (60)
and we then shall use the approximation

z(t + ∆t) + z(t -∆t) 2z(∆t). (61)
Thus, we approximate the right side of (59) by solving (60) for 0 ≤ τ ≤ ∆t, and then use (61) to compute z(t + ∆t). The last term on the right of (60) explicitly depends on z(τ ), but it involves only the unknowns associated with the nonzero entries in P . If those nonzero entries occupy only a small fraction of all unknowns, the additional effort from solving (60) will be small. Clearly, in doing so we must ensure that the overall numerical scheme remains 2sth-order accurate in time. In summary, the local time-stepping algorithm of order 2s for the solution of ( 23) computes z n+1 z(t + ∆t), given z n and z n-1 , as follows: Steps 1-3 compute the numerical solution of (60) at time τ = ∆t using the 2s-th order modified equation scheme with local time-step ∆τ = ∆t/p. Note that we recover the standard modified equation scheme by setting P = 0. Clearly, this algorithm only requires s multiplications by A(I -P ) and ps multiplications by AP . No multiplication by A is needed to compute A m-1 z n in Step 1, since

A m-1 z n = A(I -P )A m-1 z n + AP A m-1 z n = w m + u m .
We now rewrite the above alogrithm in "leap-frog manner" to determine the accuracy and establish the stability of the above algorithm. Proposition 11. The local time-stepping method (Algorithm 10) of order 2s is equivalent to

z n+1 = 2z n -z n-1 -∆t 2 A p z n , (62)
where A p is defined by

A p = A -2 s j=2 ∆t 2(j-1) (2j)! A j - 2 p 2 (p-1)s j=1 ∆t p 2(s-1+j) γ p j (AP ) j A s .
This scheme is 2s-th order accurate. Furthermore the matrix A s-1 A p is symmetric.

Proof. We do not detail the proof of the first part of the proposition, which is similar to that of Prop. 4 and 8. To prove consistency, we rewrite the scheme as

z n+1 -2z n -z n-1 ∆t 2 -2 s j=2 ∆t 2(j-1) (2j)! A j - 2 p 2 (p-1)s j=1 ∆t p 2(s-1+j) γ p j (AP ) j A s = -Az n . Since z(t + ∆t) -2z(t) + z(t -∆t) ∆t 2 -2 s j=2 ∆t 2(j-1) (2j)! A j = z (t) + O(∆t 2s )
and

(p-1)s j=1 ∆t p 2(s-1+j) γ p j (AP ) j A s = O(∆t 2s ),
Algorithm 10 is a 2s-th order approximation of [START_REF] Gustafsson | Time compact high order difference methods for wave propagation[END_REF].

As the matrices (AP ) j A s are generally not symmetric, the matrix A p is not symmetric either. However, since (A s-1 (AP ) j A s ) = A s (P A) j A s-1 = A s-1 (AP ) j A s , j ≥ 1, the matrix A s-1 A p is indeed symmetric . the two outer intervals correspond to the coarse region and the inner interval, [2 ; 4], to the refined region -refer to Section 3.3 for further details.

For every time-step ∆t, we shall take p ≥ 2 local steps of size ∆τ = ∆t/p in the refined region, with the second-order local time-stepping Algorithm 2. We choose an overlap of two and set in all instances ∆t = ∆t LF , the largest time-step allowed by the leap-frog method on an equidistant mesh with mesh size h = h coarse .

First, we consider an IP-DG discretization with P 1 elements and (small) penalty parameter α = 2, as described in Section 3.3, and the sequence of meshes T h , h coarse = 0.1, 0.05, 0.025, 0.0125, 0.00625. As we systematically reduce the global mesh size, h coarse , while simultaneously reducing ∆t, we monitor the L 2 space-time error in the numerical solution, u -u ex L 2 (0,T ;L 2 (Ω)) until the final time T = 60. In the left frame in Fig. 5, the numerical error is shown vs. the mesh size, h = h coarse . Regardless of the number of local time-steps p = 2, 4, 8, 13, the numerical method converges with order two.

We now repeat the same experiment with P 1 continuous finite elements with mass lumping for the same sequence of meshes. As shown in the right frame of Fig. 5, the local time-stepping method again yields overall second-order convergence independently of p.

Next, we consider the fourth-order time-stepping scheme (Algorithm 6) and combine it either with a continuous FE or the IP-DG discretization with P 3 elements. Thus, we expect both numerical schemes to exhibit overall fourth-order convergence with respect to the L 2 norm. Again, we choose an overlap of two and let ∆t = ∆t M E , the largest possible timestep allowed by the modified equation approach on an equidistant mesh with h = h coarse . In Fig. 6 we display the space-time L 2 -errors of the numerical solutions for the sequence of meshes h coarse = 0.2, 0.1, 0.05, 0.025, 0.0125 and different values of p. Both the continuous FE method with mass lumping and the IP-DG method, here with α = 7, yield the expected fourth-order convergence.

Finally, to validate the order of convergence of the sixth order time-stepping scheme (Algorithm 10 with s = 3), we consider the IP-DG method with P 5 elements, where we set α = 16. As above we choose an overlap of two elements and set the time-step to its (maximal) optimal value, ∆t = ∆t M E , the largest possible time-step allowed by the modified equation approach of order six on an equidistant mesh with h = h coarse . Again the numerical results shown in Fig. 7 for p = 2 corroborate the expected sixth order of convergence.

Two-dimensional example

To illustrate the usefulness of the local time-stepping methods presented above, we consider the wave equation ( 1)-( 4) with constant speed c = 1 in a computational domain Ω, that consists of two 1 × 0.95 rectangles connected by a very narrow 0.004 × 0.1 channel -see Fig. 8. We impose homogeneous Neumann conditions on the boundary of Ω and choose as initial conditions

u 0 (x) = exp( x -x 0 2 /r 2 ), x -x 0 ≤ √ 2 r, 0 otherwise (67) v 0 (x) = 0, (68)
where x 0 = (0, 0.25) and r = 0.025.

For the spatial discretization we opt for the IP-DG method with P 3 triangular elements and α = 11. Hence to resolve the propagating wave with sufficient accuracy, we find that a mesh size h coarse = 0.0125 is approximately needed. However, such triangles do not even fit inside the narrow gap, which requires h fine h coarse /16.4 to resolve its geometric features, as shown in Fig. 9.

For the time discretization we choose the fourth-order local time-stepping method from Section 4.1. Thus, the numerical method is fourth-order accurate in both space and time under the CFL restriction ∆t ≤ 0.14h, determined experimentally. If the same (global) timestep ∆t were used everywhere in Ω, it would need to be about seventeen times smaller than necessary in most of Ω, for stability reasons only. Instead, we shall use the fourth-order local time-stepping with p = 17, which for every time step ∆t = 0.14h coarse takes seventeen local time steps ∆τ = ∆t/17 inside the highly refined region.

In dimension two, or higher, the boundary between the fine and coarse mesh is typically not as well-defined as in one space dimension while the transition between larger and smaller elements is more gradual. Here to define the fine mesh, that is the subregion of Ω where local time steps are required, we surround the narrow channel by a small [-0.042 ; 0.042] × [-0.095 ; 0.095] box, as shown in Fig. 8. Those triangles, whose center of gravity lies inside that box, belong to the fine mesh; the corresponding degrees of freedom in the finite element solution are then selected merely by setting the corresponding diagonal entries of the matrix P to one -see Section 3.1.

Figure 9: The triangular mesh at various magnification rates: the darker triangles belong to the "fine mesh" and the lighter triangles to the "coarse mesh".

In Fig. 10 snapshots of the numerical solution u are shown at times t = 0.02, 0.67, 0.13, 0.17, 0.24, 0.27. A circular wave is initiated by the Gaussian pulse centered about x 0 in the upper region, which propagates outward until it impinges on the lower boundary at t = 0.67.

Then, a fraction of the wave penetrates the channel and generates a circular outgoing wave as it reaches the opposite lower region. Further reflections occur as the wave moves back and forth inside the channel, subsequently generating multiple circular waves in the upper and lower domains. Finally, we verify that the discrete numerical energy defined in (56) is truly conserved over time. In Fig. 11 we follow the time evolution of the energy and its relative variation, which remains within machine precision for the entire simulation until T = 8.7, that is during 5000 time-steps.

Conclusion

We have presented explicit local time-stepping methods for the wave equation, which allow arbitrarily small time-steps precisely where the smallest elements in the mesh are located. When combined with a symmetric finite element discretization in space with an essentially diagonal mass matrix, the resulting discrete time-marching scheme remains truly explicit, while it also conserves a discrete energy. Starting from the standard second-order "leap-frog" scheme, we have derived local time integration methods of arbitrary order. When the "fine" region, where local time-steps are used, slightly extends into the surrounding "coarse" region of the mesh, we find that the resulting numerical scheme permits the use of the optimal maximal time-step, dictated by the coarse mesh size.

Since the local time-stepping methods presented here are truly explicit, their parallel implementation is straightforward. Let ∆t denote the time-step imposed by the CFL condition in the coarser part of the mesh. Then, during every (global) time-step ∆t, each local time-step of size ∆t/p inside the fine region of the mesh, with p ≥ 2 any integer, simply corresponds to sparse matrix-vector multiplications that only involve the degrees of freedom associated with the fine region of the mesh. Those "fine" degrees of freedom can be selected individually and without any restriction by setting the corresponding entries in the diagonal projection matrix P to one; in particular, no adjacency or coherence in the numbering of the degrees of freedom is assumed. Hence the implementation is straightforward and requires no special data structures. The local time-stepping methods derived here for the scalar wave equation immediately apply to more general second-order hyperbolic problems, as in elasticity or electromagnetics, for which either symmetric discontinuous Galerkin [START_REF] Grote | Interior penalty discontinuous Galerkin method for Maxwell's equations: energy norm error estimates[END_REF][START_REF] Grote | Interior penalty discontinuous Galerkin method for Maxwell's equations: optimal L 2 -norm error estimates[END_REF] or mass lumping techniques are available [START_REF] Elmkies | Finite elements and mass lumping for Maxwell's equations: the 2D case[END_REF]. They also generalize to the situation of nonzero forcing. Clearly these timestepping schemes can also be combined with finite difference methods on highly stretched grids, such as body-fitted grids [START_REF] Villamizar | Time-dependent numerical method with boundary-conforming curvilinear coordinates applied to wave interactions with prototypical antennas[END_REF], if the underlying finite difference discretization leads to a symmetric stiffness matrix. In the presence of hierarchical mesh refinement, each local timestep in the fine region can itself include further local time-steps inside a smaller subregion with an even higher degree of local mesh refinement.

Algorithm 2 . 1 . 2 ∆t p 2 ( 2 w

 21222 Set w = A(I -P )z n and z0 = z n ; 2. Compute z1/p = z0 -1 w + AP z0 ); 3. For m = 1, . . . , p -1, compute z(m+1)/p = 2z m/p -z(m-1)/p -∆t p + AP zm/p ; (36) 4. Compute z n+1 = -z n-1 + 2z 1 .

Figure 3 :

 3 Figure3: IP-DG P 1 elements. The eigenvalues of (∆t 2 /4)A p : overlap by one element (left); overlap by two elements (right). The vertical scale is strongly magnified: 0.99999 < λ max < 1.00001. elements (IP-DG P 1 elements).

Algorithm 10 . 1 . 3 .(- 1 ) 4 .

 101314 Set z0 = z n , w m = A(I -P )A m-1 z n , u m = AP A m-1 z n , m = 1, . . . , s; 2. Compute v m = (-1) m (w m + u m z0 ), m = 1, . . . , s; For m = 1, . . . , p -1, compute v 1 = -w 1 + s l=1 (-1) l (2l)! w l -AP zm/p , k = 1, . . . , s; v k = (-1) k w k + s l=k+1 l-k (2(l -k))! w l -AP v k-1 , k = 2, . . . , s;z(m+1)/p = 2z m/p -z(m-1)/p + Compute z n+1 = -z n-1 + 2z 1 .

Figure 5 :

 5 Figure 5: Error vs. h = h coarse for P 1 finite elements with p = 2, 4, 8 and 13: IP-DG (left) and continuous FE (right).

Figure 6 :

 6 Figure 6: Error vs. h = h coarse for P 3 finite elements with p = 2, 4, 6, 7: IP-DG (left) and continuous FE (right).

Figure 7 :

 7 Figure 7: Error vs. h = h coarse for the IP-DG method with P 5 elements with p = 2.

Figure 8 :

 8 Figure 8: Two-dimensional example: the computational domain Ω.

Figure 10 :

 10 Figure 10: Two-dimensional example: the solution is shown at times t=0.02, 0.67, 0.13, 0.17, 0.24, 0.27.

Figure 11 :

 11 Figure 11: Time evolution of the discrete energy defined in (56) (left) and of its relative variation |E(t)/E(0) -1| (right).

  1 2 in (43), which guarantees stability if and only if E n+ 1 2 is positive semidefinite or, equivalently, if and only if the matrices I -∆t 2 4 A p and A p are both positive semidefinite. Hence if λ min and λ max denote the smallest and largest eigenvalues of A p , respectively, the numerical scheme will be stable if and only if

	0 ≤	∆t 2 4	λ min ≤	∆t 2 4	λ max ≤ 1.

Table 1 :

 1 IP-DG P 1 elements. The largest eigenvalue of (∆t 2 /4)A p with ∆t = ∆t opt for an overlap by one element.

	p

  IP-DG with P 1 elements. The eigenvalues of (∆t 2 /4)A p : without overlap (left); with an overlap by one element (right).

	1.6									1.6
	1.4									1.4
	1.2									1.2
	1									1
	0.8									0.8
	0.6									0.6
	0.4									0.4
	0.2									0.2
	0									0
	0 -0.2	0.2	0.4	∆t/∆t opt	0.6	0.8	1	0 -0.2	0.2	0.4	∆t/∆t	0.6	0.8	1
	1.00001									1.00001
	1									1
	0.7 0.99999	0.75	0.8	∆t/∆t opt	0.85	0.9	0.95	1	0.7 0.99999	0.75	0.8	∆t/∆t opt	0.85	0.9	0.95	1

opt Figure 2:

Table 2 :

 2 IP-DG P 1 elements. The largest eigenvalue of (∆t 2 /4)A p with ∆t = ∆t LF for an overlap by two elements.

	p

Table 3 :

 3 Continuous P 1 -FE. The largest eigenvalue of ∆t 2 4 A p at ∆t = ∆t LF for an overlap by

	one element.
				p
	h coarse	2	3	4	10	13
	0.5	0.9828 0.9792 0.9993 0.9999 0.9999
	0.2	0.9969 0.9962 0.9999 0.9999 0.9999
	0.1	0.9992 0.9991 0.9999 0.9999 0.9999
	0.05	0.9998 0.9998 0.9999 0.9999 0.9999
	0.025 0.9999 0.9999 0.9999 0.9999 0.9999

Table 4 :

 4 IP-DG P 3 elements. The smallest eigenvalue of AA p and A -(∆t 2 /4)AA p at ∆t = ∆t M E with p = 2: without overlap (left); with overlap by one element (right).

	h coarse 0.5 0.2 0.1 0.05 0.025 -7.10 10 AA p -5.10 5 -2.10 7 -3.10 8 -5.10 9	A -(∆t 2 /4)A p 1.1 1.1 1.1 1.1 1.1

Table 5 :

 5 P 3 FE elements. The smallest eigenvalue of AA p and A-(∆t 2 /4)AA p for ∆t = ∆t M E with p = 2 without overlap. h coarse AA p A -(∆t 2 /4)AA p

	0.5	1.2	0.15
	0.2	1.2	0.9
	0.1	1.2	0.53
	0.05	1.2	1.1
	0.025	1.2	1.1

RR n 6377

Acknowledgements We thank Ernst Hairer and Christian Lubich for useful comments and suggestions.

Energy conservation and stability

We shall now exhibit a discrete energy which is conserved by the fourth-order local timestepping scheme (Algorithm 6) and thereby determine a necessary and sufficient condition for stability. Next, we perform a systematic numerical study to demonstrate the improvement in the CFL condition achieved by slightly extending the overlap of the fine into the coarse region.

First, using Proposition 8 we rewrite the local time-stepping scheme as (54)-(55). As a consequence, we can prove the conservation of the following energy.

Proposition 9. The fourth-order local time-stepping scheme (Algorithm ( 6)) conserves the following energy:

Proof. As A p is not symmetric, we first have to pre-multiply (54) by A. Next, we take the inner product of the resulting expression with z n+1 -z n-1 , which yields the conservation of

That expression is equal to E n+ 1 2 above, since

Since both A and AA p are symmetric matrices, E n+ 1 2 is a symmetric quadratic form which is positive for ∆t sufficiently small. Hence the local time-stepping scheme will be stable if, and only if, the two matrices A -∆t 2 4 AA p and AA p are positive semidefinite. For p = 1, we have AA p = A A -∆t 12 A 2 and thus we recover the CFL condition of the modified equation scheme:

where λ max is the largest eigenvalue of A.

For p ≥ 2, the CFL condition cannot be determined analytically. Thus, we shall again provide a systematic numerical study to evaluate the stability of Algorithm 6 in a typical situation.

We consider the same one-dimensional problem as in Section 3.3, but now use an IP-DG discretization with P 3 -elements with (small) penalization α = 7 to achieve fourth-order accuracy both in space and time. In the absence of local time-stepping, that is with p = 1, the maximal time-step dictated by the equidistant (coarse) mesh is ∆t M E = 0.265 h coarse . For p ≥ 2, we choose the same optimal time-step, ∆t = ∆t M E and compute the smallest eigenvalue of A -∆t 2 4 AA p and AA p , respectively -both need to be positive for stability.

We can now prove the conservation of a discrete energy.

Proposition 12. The local time-stepping scheme (Algorithm 10) of order 2s conserves the following energy:

Proof. As A p is not symmetric, we first pre-multiply (62) by A s-1 . Next we take the inner product of the resulting expression by z n+1 -z n-1 , which yields the conservation of

That expression is equal to E n+ 1 2 above, since

Since both A s-1 and A s-1 A p are symmetric matrices, E n+ 1 2 is a symmetric quadratic form which is positive for ∆t sufficiently small. Hence the local time-stepping scheme will be stable if, and only if, the two matrices A s-1 -∆t 2 4 A s-1 A p and A s-1 A p are positive definite.

Numerical Results

We shall now present numerical experiments that confirm the expected order of convergence and demonstrate the versatility of the above local time-stepping methods. First, we consider a simple one-dimensional test problem to show that the different local time-stepping schemes presented above indeed yield the expected overall rate of convergence when combined with a spatial finite element discretization of comparable accuracy, independently of the number of local time-steps p. Then, we consider wave propagation in two space dimensions with a locally highly refined mesh to illustrate the usefulness of local time-stepping in the presence of complex geometry.

Convergence study

We consider the one-dimensional wave equation with constant wave speed c = 1 on the interval Ω = [0 ; 6] with periodic boundary conditions. The initial conditions are chosen to yield the exact solution u ex (x, t) = sin(8π(x -t)/3), which corresponds to a sinusoidal wave propagating to the right. Again, we divide Ω into three equal parts. The left and right intervals, [0 ; 2] and [4 ; 6], are discretized with an equidistant mesh of size h coarse , whereas the interval [2 ; 4] is discretized with an equidistant mesh of size h fine = h coarse /p. Hence, INRIA


Pascale Henner

Michel Schiavon

Jean-Louis Morel

Eric Lichtfouse
email: eric.lichtfouse@dijon.inra.fr

Polycyclic aromatic hydrocarbons (PAH) occurrence and remediation methods

Polycyclic aromatic hydrocarbons (PAHs) are potentially mutagenic and carcinogenic substances occurring at various concentrations in atmosphere, soils, waters and sediments. PAHs, inherited both from natural and anthropogenic processes, are persistent organic pollutants (POP) due to their chemical stability and biodegradation resistance. The increase of road transportation, and of industrial and agricultural activities has led to a notable build up of PAH amounts in the environmental media. For cases of heavy soil pollution, the main remediation methods are containment, thermal desorption, incineration, and microbial degradation. Phytoremediation, a low-cost process based on microbial activation in the root zone, is a novel method under deep investigation.

Résumé

Les hydrocarbures aromatiques polycycliques (HAP) sont des substances potentiellement mutagènes et cancérigènes présentes à des concentrations variables dans l'atmosphère, les sols, les eaux et les sédiments. Les HAP, issus de processus naturels ou anthropiques, sont des polluants organiques persistants (POP) de part leur inertie chimique et leur résistance à la biodégradation. L'accroissement des transports routiers, des activités industrielles et agricoles a induit une augmentation notable de la quantité d'HAP dans l'environnement moderne. Dans les situations de pollutions massives du sol, les principales techniques de réhabilitation sont l'isolement, la désorption thermique, l'incinération, et la biodégradation microbienne. La phytoremédiation des HAP, un processus bon marché fondé sur l'activation des micro-organismes du sol par les racines, est une nouvelle méthode en plein développement.

Introduction

Polycyclic aromatic hydrocarbons (PAHs), some of which are known precursors of mutagenic derivatives, are widely occurring in natural media such as soil, sediment, water, air and plants as a result of both natural and anthropogenic processes [1 and refs. therein]. The main PAH sources are highlighted in Table 1. Since natural fires and petroleum formation have been occurring all through earth times, PAHs have undoubtedly been circulating through biogeochemical cycles for millions of years. Nonetheless, the increasing use of fossil fuels and their combustion products by human beings during the two past centuries raises several questions about PAHs hazards for living organisms. First, apart from accidental oil spills leading to massive pollutions, the precise origin of trace PAHs, e.g. natural vs. anthropogenic, is rarely clear. Second, the toxicity of PAHs, like other hazardous chemicals, requires their bioavailability. And since most PAHs are highly hydrophobic, their pathways of transfer through geological and biological media are far from being understood. Third, unambigous correlations between PAH sources and carcinogenic effects have been reported only for intense exposure to PAHs such as for coal-mine workers. Here we report the major aspects about PAHs in the literature with emphasis on molecular structures and properties, origin and distribution, analytical methods, fate in the soil-plant system and remediation techniques.

Molecular structure and properties

The molecular structures of the 16 PAHs defined as priority pollutants by the American Environmental Protection Agency (EPA) are drawn on Figure 1. PAHs are relatively neutral and stable molecules, including one or more benzenic rings, substituted or not, on a linear or angular layout. PAHs have low solubilities [2,3] and low volatilities, except small components like naphthalene. Their lipophilicity is high, as measured by water-octanol partition coefficients (Kow). PAH show long half-lives in geological media. In an aerobic sediment, for example, half lives range from 3 weeks for naphthalene up to 300 weeks for benzo[a]pyrene [4]. PAHs are regarded as persistent organic pollutants (POP) in the environment. This persistence is increasing with ring number and condensation degree [5].

Toxicity

The carcinogenicity levels of some PAHs are given in Table 2. Monoaromatic hydrocarbons such as benzene-toluene-xylene mixtures (BTX) are known phytotoxic [6], albeit the mechanisms of growth inhibition are unclear. BTX either disturb exchanges between soil and atmosphere, or they alter plant tissues . Higher PAHs are less phytotoxic in the short term, but some of them and their metabolites are suspected to be mutagenic and carcinogenic (Table 1). Indeed, their ability to form DNA adducts can induce mutations.

Table 2. Carcinogenicity levels of PAHs [7]. Keys : (-) no effect , (+-) very weakly active, (+) weakly active, (+++) active, (++++) very active.

________________________________________ PAH CARCINOGENICITY ________________________________________ Naphthalene (-) Anthracene (-) Pyrene (-) Chrysene (+-) Benz[a]anthracene (+) Dibenz[a,h]anthracene (+++) Benzo[a]pyrene
(++++) ________________________________________ Some PAHs, like benzo[a]pyrene, are immunodepressive in a dose-response manner, a property correlated with carcinogenicity [7]. PAHs present a risk for the environment, especially if they can enter the food chain through contact, inhalation or ingestion. A link has been disclosed between coal exposure and cancer among coal miners and chimney sweepers [7], but here other carcinogenic substances than PAHs may be responsible. In fact, very few data assessing a correlation between a pure PAH and cancer in humans are available. Effects on human beings are therefore assumed to be the same as the ones observed on mice, guinea-pigs, and other animals.

Indigenous sources

PAHs may be indigenous as a result of rapid aromatization of biological compounds during humification processes. They can also reflect contributions from ancient hydrocarbon sources, such as coal residues, provided that the soil has been initially formed by weathering of a petroleum source-rock or coal seam. These sources are regarded as natural or pristine.

Exogenous sources

These sources include either natural inputs, e.g. the migration of petroleum from deeper horizons, or anthropogenic inputs such as fossil fuel spills and deposition of airborne particles from car exhausts. Nevertheless, the precise origin of PAHs in airborne particles from pristine areas is difficult to assess. Vegetation fires are both natural and anthropogenic sources of PAHs. It is sometimes possible to distinguish, on structural grounds, petrogenic sources that account for the direct contribution of uncombusted fossil fuel, from pyrolytic or pyrogenic sources formed at relatively high temperatures by the incomplete combustion of fossil fuels and vegetation fires.

Distribution

PAHs are found worldwide in soils and sediments, even in the most pristine areas such as the Poles [8][9][START_REF] Tarradellas | Soil Ecotoxicology[END_REF]. Reported background contents of PAHs in soils of North America and Europe extend from 50 to 500 µg.kg -1 [9]. Because vehicles, heating systems and waste incineration exhausts are the main anthropogenic sources, PAHs contents increase in the vicinity of large urban and industrial areas [4,8,[START_REF] Tarradellas | Soil Ecotoxicology[END_REF]. PAHs are carried at long distances by air in the troposphere, bound to particulate carbon. Thus, they can form a worldwide deposit by wet or dry deposition [START_REF] Youngblood | [END_REF]. Indeed, the qualitative distribution of PAHs in modern media are often similar, strengthening the influence of the well-mixed airborne source [8]. This source include a major part of PAHs produced by pyrolysis at moderate temperature, e.g. fires and car exhausts. Since the beginning of the industrial revolution, the rate of deposition has exceeded the rate of degradation in the soils [4]. Consequently, PAHs accumulate in the top soil.

Fate of PAHs

The main pathways of PAHs are shown on Figure 2. PAHs are hydrophobic compounds whose mobility in the soil-plant system is low [3]. Equilibriums between solid, aqueous and vapour phases are very slow. As for pesticides [12], PAHs and their metabolites are expelled from the aqueous phase to adsorb on the hydrophobic surfaces such as organic matter [13][14][15]. They are either trapped in the pores, fixed with covalent or hydrogen bonds, or bound during humification processes [16]. PAHs structure and stability stand in the way of their biodegradation [17] by micro-organisms, e.g. fungi and bacteria. Biodegradation is slow and is a function of environmental parameters such as oxygen, water and nutriment contents. Migration of PAHs from the top soil is slow. PAHs seem to migrate bounded to particles [4]. The major ways of entry of PAHs in the plant seems to be through the leaves, from the vapour phase and by contact with contaminated soil particles [18]. Entry through roots and translocation to shoots is still hypothetic because the way to xylem requires PAHs (compounds with very low solubilities) to solubilize in the symplasm. PAHs seem also to adsorb on the root cell walls.

Analysis

Extraction of PAHs in soils is mainly conducted with several methods : shake, humic acid extraction, supercritical fluid, soxhlet and ultrasonic bath. Soxhlet and ultrasonic extractions give similar yields [19]. Yields decrease with increasing soil organic matter and polluant contents. Bound residues are usually not available with these techniques [START_REF] Mahro | In : Bioremediation of chlorinated and polycyclic aromatic hydrocarbon compounds[END_REF][START_REF] Kästner | [END_REF]. A part of bound residues may be recovered by humic acid extraction, e.g. with 0.1 M NaOH. Identification and quantification of PAHs is currently achieved by liquid chromatography coupled to UV detector, or by gas chromatography coupled to mass spectrometry (GC-MS).

Remediation of PAH contaminated soils

The strategies of PAH remediation depend on the extent of contamination, the horizontal and vertical limits of the pollution as well as its qualitative and quantitative content [22]. While containment withdraws temporarily the pollutants, the other treatments attempt to remove, stabilize or degrade them.

Containment

Contaminated materials are enclosed in physical, chemical or hydraulic barriers which prevent the polluants to migrate and contamine groundwaters [START_REF]ADEME Les différentes techniques de réhabilitation des sols et sites pollués[END_REF][START_REF] Cunningham | In Vitro Cell[END_REF][25]. In France, monitoring rules for containment are as strong as for rubbish dumps. For example, waters circulating through the site must be continuously collected and analysed. Access to and use of the site are restricted. This technique is usually temporarily used while waiting for an another way of treatment.

Thermal techniques

Incineration and thermal desorption are the most effective and used techniques, although their cost is relatively high because they involve soil excavation, transport then heat treatment [22][START_REF]ADEME Les différentes techniques de réhabilitation des sols et sites pollués[END_REF][START_REF] Cunningham | In Vitro Cell[END_REF][25]. Soil can be used freely after treatment.

Physical and chemical techniques

Polluants may be first extracted then transfered in gaz or liquid phases. These fluids are then specifically treated. For instance, volatile and soluble compounds can be extracted with vacuum, air or water. Non volatile or non soluble compounds such as high molecular weight PAHs are extracted with the help of surfactants, acids, heat or physical treatments to enhance their availibility [22][START_REF]ADEME Les différentes techniques de réhabilitation des sols et sites pollués[END_REF][START_REF] Cunningham | In Vitro Cell[END_REF][25]. Alternatively, when the later methods fail, the application of an electric current may favor the migration of pollutants by electro-osmosis.

Microbial degradation

Bioremediation uses the natural capacity of soil micro-organisms to degrade PAHs. Degradation conditions are optimized by aeration, agitation, moistening and addition of nutriments. Genetically engineered strains, with multiple degradative abilities, are usually less competitive that well-adapted natural strains. Landfarming is a technique whereby polluants are incorporated in a top layer of an agricultural soil [5,26]. No culture is conducted on the soil during the treatment. Nutriments are incorporated to stimulate bioactivity. Degradation is performed during the whole year, except in winter when humification process is greater. Composting is another way of treatment. Composting of petroleum by-products may be performed with high contents of polluant. Here, the compost structure is optimized by adding agents like sawdust, wood cuttings, coco fibers, moss, sand or vermiculite. Enhanced degradation is achieved by aeration, humidification and regular mixing of compost piles. On the other hand, bioreactors are laboratory scale units now effectively used on an industrial scale. Bioreactors are either fixed or rotative. They can be loaded either with soil samples or with slurry phases [START_REF] Mahro | In : Bioremediation of chlorinated and polycyclic aromatic hydrocarbon compounds[END_REF][27][28][29][30]. Here, the inoculum is usually either a given volume of soil with its natural microflora or engineered strains. Because the degradative ability depends on plasmids, new plasmid-enriched strains have been designed.

Phytoremediation

Phytoremediation is the use of plants, their associated microflora and agronomical techniques to lower soil toxicity by extraction, degradation and stabilisation of the pollutants [22]. It fullfils the need of an efficient and low-cost on-site remediation technique. Here, plant roots modify favourably the life habitat of micro-organisms by exudate release [31,32], nutrients supply, e.g. nitrogen and sugars, and aeration and water drainage [22]. Roots also increase the contact surface between micro-organisms, soil and pollutants [START_REF] Qiu | In : Bioremediation through rhizosphere technology[END_REF]. Further, the number and diversity of microorganisms are greatly increased in the rhizosphere, as well as the metabolic capacities and specific associations like myccorhizae [START_REF] Qiu | In : Bioremediation through rhizosphere technology[END_REF]. Therefore, biodegradation of PAHs should be quicker and greater in the presence of plants.

Indeed, several recent studies show that the concentration of organic pollutants such as PAHs can be notably reduced by vegetation [13,31,[START_REF] Qiu | In : Bioremediation through rhizosphere technology[END_REF][34][35][START_REF] Schwab | Bioremediation through rhizosphere technology[END_REF][37]. However, since these results are based on the observed decrease of the parent molecule, the molecular mechanisms of dissipation are still unclear. It is not known, for example, whether plants could drain PAHs and their metabolites upward to shoots and leaves. Dissipation may also results from enhanced formation of bound residues which could be released in the long term. Lastly, production, solubilization and transfer of toxic metabolites has to be studied in order to predict a possible migration of the pollution toward other media such as groundwater and plants. Nonetheless, phytoremediation appears as a mild and low cost method to remediate contaminated soils from ancient industries, railways, etc.

Figure 1 .

 1 Figure 1. Molecular structure of the 16 polycyclic aromatic hydrocarbons (PAHs) selected as priority pollutants by the American Environmental Protection Agency (EPA). Numbers are molecular weights.

Figure 2 .

 2 Figure 2. Main possible pathways of PAHs in the soil-plant system.

Table 1

 1

* Author for correspondance