Frédérique Le Louër

A high order spectral algorithm for elastic obstacle scattering in three dimensions

Keywords: Navier equations, Dirichlet condition, Neumann condition, boundary integral equations, spectral method

HAL is

Introduction

In this paper we present a high order spectral algorithm for solving the time-harmonic Navier equations in the exterior of a three dimensional bounded obstacle, with a Dirichlet or a Neumann boundary condition. Efficient numerical solution of such scattering problems is of practical interest for various industrial applications related to inverse scattering problems such as non-destructive testing of materials [START_REF] Bonnet | Inverse problems in elasticity[END_REF].

The method of boundary integral equation (BIE) is a classical tool for solving scattering problems of timeharmonic waves in unbounded, homogeneous and isotropic media. Its efficiency has been widely demonstrated for low and middle frequency waves via boundary collocation methods or boundary element methods (BEM). However, the boundary integral formulation has the drawback that it reduces the time-harmonic scattering problem to a fully-populated and, in general, non-Hermitian complex linear system. The complexity of direct solvers to invert the scattering matrix is O(N 3 ) in time and O(N 2 ) in memory if N is the number of degrees of freedom. The complexity of iterative solvers such as GMRES, which iteratively build a sequence of approximate solutions, is O(n it × N 2 ) in time and memory if n it is the number of iteration. The boundary element mesh size depends on the frequency of the problem and the geometry of the obstacle. Therefore, the solution of realistic problems is limited by the number of degrees of freedom that can be handled by one or the other solver on any machine. In the scientific literature, there exist various methods to decrease the computation time and the memory size: preconditionning and fast methods. Preconditionning methods aim at reducing significantly the number of iteration of the iterative solver applied to a full system-matrix. Fast methods aim at reducing the size of the discretization space. One of them is the so-called fast multipole method (FMM) which speed up the matrix-vector multiplication used by iterative solvers with BEM discretization. Earlier works in acoustics and electronagnetism lead to a complexity of O(N log 2 N ) in time and memory per iteration. It has recently been extended to linear elasticity by Bonnet, Chaillat and Semblat [START_REF] Chaillat | A fast multipole method formulation for 3d elastodynamics in the frequency domain[END_REF][START_REF] Chaillat | A multi-level fast multipole bem for 3-d elastodynamics in the frequency domain[END_REF][START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF]. A different kind of compression for the BE matrices can be obtained by applying the adaptive cross approximation (ACA) algorithm to hierarchical matrices [START_REF] Bebendorf | Approximation of boundary element matrices[END_REF][START_REF] Bebendorf | Accelerating Galerkin BEM for linear elasticity using adaptive cross approximation[END_REF][START_REF] Bebendorf | Adaptive low-rank approximation of collocation matrices[END_REF]. A comparison between these two techniques can be found in [START_REF] Brunner | Comparison of the fast multipole method with hierarchical matrices for the Helmholtz-BEM[END_REF] in the low-and high-frequency regimes for the Helmholtz equation. Other fast methods, as spectral methods, accelerate the convergence of the discrete solution so that the BIE can be reduced to linear systems with substantially fewer unknowns than with other BEM-based methods.

It is the purpose of this paper to use high order spectral methods to solve the BIEs in linear elasticity, posed on simply connected closed surfaces. It was first introduced by Ganesh and Graham for acoustic scattering [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF] and consists in transforming the integral equation on the surface Γ to an integral equation on the unit sphere using a change of variable and then by expanding the integrand and looking for a scalar solution in terms of series of scalar spherical harmonics. The integral equation is finally reduced to a linear system whose matrix coefficients are the integral equation operator evaluated at each basis functions. Superalgebraic convergence of the discrete solution was proven for smooth obstacles and Fredholm integral equations of the second kind.

This technique has then been extended to electromagnetism by Ganesh and Hawkins for solving the perfect conductor problem [START_REF] Ganesh | A spectrally accurate algorithm for electromagnetic scattering in three dimensions[END_REF] by expanding the integrand and looking for a solution in terms of (componentwise) series of scalar spherical harmonics. Spectral accuracy of the algorithm was proven for the magnetic field integral equation (MFIE) only. In [START_REF]An efficient surface integral equation method for the time-harmonic Maxwell equations[END_REF] they exploited the vector spherical harmonics, but this approach does not yield spectrally accurate approximations : a stagnated convergence is observed (see [START_REF] Ganesh | A spectrally accurate algorithm for electromagnetic scattering in three dimensions[END_REF]Remark 2]). However, as being eigenfunctions of integral operators on the sphere, vector spherical harmonics are a natural choice for surface formulation of three dimensional scattering problems. Expanding the integrand in terms of componentwise sum of scalar spherical harmonics, one obtains the hybrid high-order method [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF] with a complexity of O(n 5 ), where n is the largest degree of the spherical harmonics. For smooth obstacles, the spectral convergence rate is O(n -s ), for any positive integer s. The order constant in the spectral accuracy depends on the size, the geometry, and the smoothness of the obstacles.

Further efficient algorithms for non-spherical perfect conductors require construction of high-order basis functions that are also tangential on the surface of the non-spherical scatterers. This was achieved by Ganesh and Hawkins in [START_REF]A high-order tangential basis algorithm for electromagnetic scattering by curved surfaces[END_REF], concluding the development of high-order algorithm for computation of tangential surface currents for electromagnetic scattering by a single perfect conductor. For elastic obstacle scattering the unknown field is not required to be tangential and hence in this work, we follow the approach developed in [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF].

It is well known that BIEs in linear elasticity involve integral operators with a strongly singular or hypersingular kernel. To apply the hybrid spectral method of Ganesh and Hawkins we use a standard technique coming from the finite element theory, that is integration by parts. Indeed, the boundary integral operators with a strongly or hypersingular kernel -that is the double layer boundary integral operator and the Neumann trace of the single and double layer potentials -can be expressed in terms of integral operators with a weakly singular kernel and surface differential operators [START_REF] Hsiao | Boundary integral equations[END_REF]. Surface derivatives of the vector spherical basis functions can be computed analytically, thus no more approximations is required. By this way, one expects that the spectral accuracy of the hybrid method developed for the MFIE naturally extends to the BIEs of the first kind in linear elasticity.

In the low and medium frequency region, the resulting full matrix system is sufficiently small (few thousands of unknowns) to be solved with a direct solver. To speed up the implementation of the arrays, one can use the Fast Fourier Transform. For higher frequencies and complex geometries (in terms of non-convexity and angularity) an iterative solver is required. Ganesh and Hawkins developed in [START_REF]Sparse preconditioners for dense complex linear systems arising in some radar cross section computations[END_REF] an algebraic preconditioner specifically for high order spectral algorithms using sparse approximations.

Let Ω ⊂ R 3 be a bounded domain with a simply connected closed boundary Γ of class C 2 at least and outward unit normal vector n and let Ω c denote the exterior domain R 3 \Ω. Throughout the paper we denote by H s loc (Ω c ) and H s (Γ) the standard (local in the case of the exterior domain) complex valued, Hilbertian Sobolev space of order s ∈ R defined on Ω c and Γ respectively (with the convention H 0 = L 2 .) Spaces of vector functions will be denoted by boldface letters, thus

H s = (H s ) 3 . We set H 1 loc (Ω c , ∆ * ) := u ∈ H 1 loc (Ω c ) : ∆ * u ∈ L 2 loc (Ω c
) . We use the notation of [2, section 2] for the formulation of the elastic scattering problems and quote some important properties of the boundary integral operators from [15, section 2]. We assume that the Lamé parameters µ and λ + µ and the density ρ are positive constants. The propagation of time-harmonic elastic waves in the threedimensional isotropic and homogeneous elastic medium is described by the Navier equation

∆ * u + ρω 2 u = µ∆u + (λ + µ)∇ div u + ρω 2 u = 0, (1.1)
where ω > 0 is the frequency. In this work we are interested in computing an approximation to the radiating solution and its far field pattern of the following exterior boundary value problems for elastic waves: Given vector densities f ∈ H

u s = 0, r = |x|, (1.4)
uniformly in all directions. Here, the longitudinal wave is given by u p = -κ -2 p ∇ div u and the transversal wave is given by u s = uu p associated with the respective wave numbers κ p and κ s given by κ 2 p = ρω 2 (λ + 2µ) -1 and κ 2 s = ρω 2 µ -1 . The radiation condition implies that the solution has an asymptotic behavior of the form

u(x) = e iκp|x| |x| u ∞ p (x) + e iκs|x| |x| u ∞ s (x) + O 1 |x| , |x| → ∞,
uniformly in all directions x = x |x| . The fields u ∞ p and u ∞ s are defined on the unit sphere S 2 in R 3 and known as the longitudinal and the transversal far-field pattern, respectively. Setting G(κ, xy) = e iκ|x-y| 4π|x -y| , the fundamental solution of the Navier equation is given by

Φ(x, y) = 1 µ G(κ s , x -y)I R 3 + 1 κ 2 s ∇ x T ∇ x G(κ s , x -y) -G(κ p , x -y) .
It is a 3×3 matrix-valued function and we have Φ(x, y) = T Φ(x, y) = Φ(y, x). For a solution to the Navier equation (1.1) that satisfies the Kupraze radiation condition, one can derive the Somigliana integral representation formula for x ∈ Ω c :

u(x) = Γ T [P y Φ(x, y)]u(y) -Φ(x, y)P y u(y) ds(y), (1.5)
where P y = P (n(y), ∂ y ) and P y Φ(x, y) is the tensor obtained by applying the traction operator P y to each column of Φ(x, y). The elastic scattering problems can be reduced, via several different use of (1.5), to a single combined-field boundary integral equation (CFIE) [START_REF] Bonnet | Boundary integral equations methods in solids and fluids[END_REF][START_REF] Charalambopoulos | On the Fréchet differentiability of boundary integral operators in the inverse elastic scattering problem[END_REF][START_REF] Kupradze | Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF]. In Section 2 we review the standard direct and indirect approaches. We also give variant representations of the boundary integral operators that are better suitable for their numerical evaluations. These formulations are obtained with the help of the tangential Günter derivative (see [START_REF] Hsiao | Boundary integral equations[END_REF][START_REF] Kupradze | Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF] or Lemma 2.1 for a definition). Section 3 and 4 are concerned with the description of a high order spectral algorithm for solving boundary integral equations in linear elasticity. We globally use the same notations and steps as Ganesh and Hawkins in [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF]. We start in Section 3 by transporting the boundary integral operators and the surface differential operators over the unit sphere and by splitting the kernels into a weakly singular part and a smooth part. The pull-back of surface differential operators on the unit sphere are evaluated using straightforward equalities that link them to the corresponding differential operators on the unit sphere. Then in Section 4, we summarize the procedure to obtain a fully discrete approximations of the CFIEs. Numerical results in the low and medium frequency region for a variety of three-dimensional convex and non-convex smooth obstacles are presented in Section 5 to show the efficiency of the method. And especially, we compare the fast spectral algorithm (FSM) with the fast multipole method (FMM). The results demonstrate that the FSM still is very competitive with the FMM in linear elaticity as was shown by Ganesh and Graham in [17, page 234] in acoustic scattering. Finally, we draw concluding remarks and we discuss possible research lines in Section 6.

The solution of elastic obstacle scattering problems

The scattering problem of time-harmonic elastic waves by a bounded obstacle Ω leads to special cases of the above boundary value problems (1.1)-(1.4) which can be reduced in several different ways to a single combinedfield boundary integral equation. We review the standard direct and indirect approaches. For the Dirichlet problem we follow [13, Section 2] and for the Neumann problem we extend acoustic techniques [START_REF] Kress | Linear integral equations[END_REF]Chapter 5].

Let us start with the rigid body problem: the total displacement field u + u inc is given by the superposition of the incident field u inc , which we assume to be an entire solution of the Navier equation, and the scattered field u, which solves the Navier equation in Ω c , and satisfies the Dirichlet boundary condition (1.2) with f = -u inc and the Kupradze radiation condition. In this case, the direct approach is based on the Somigliana integral representation formula (1.5) and the following result for

x ∈ Ω c [15, formula (2.6)] 0 = Γ T [P y Φ(x, y)]u inc (y) -Φ(x, y)P y u inc (y) ds(y).
Adding the last equation to (1.5) and using the boundary condition of u, we obtain the following integral representation:

u(x) = - Γ Φ(x, y)P y u(y) + u inc (y) ds(y), x ∈ Ω c . ( 2
SP (u + u inc ) = 2u inc and (I + D ′ )P (u + u inc ) = 2P u inc on Γ,
where S is the single layer boundary integral operator and D ′ is the traction derivative of the single layer potential, defined by

Sϕ(x) = 2 Γ Φ(x, y)ϕ(y) ds(y), D ′ ϕ(x) = 2 Γ [P x Φ(x, y)]ϕ(y) ds(y).
Setting ϕ = P (u + u inc ), the field u given by (2.1) solves the Dirichlet boundary value problem (1.1)-(1.2) if the density ϕ solves the following integral equation

ϕ + D ′ ϕ + iηSϕ = 2 P u inc + iηu inc on Γ. (2.2)
where η is a non-zero real constant. The operator S :

H 1 2 (Γ) → H 1 2 (Γ) is compact while the operator D : H 1 2 (Γ) → H 1 2 (
Γ) has a strongly singular kernel. It can be shown that the homogenous form of (2.2) only has the solution ϕ = 0, however we cannot use immediately Riesz theory to prove existence of a solution to (2.2). We use a left equivalent regularizer [27, Chapter IV, Section 5] to modify the integral equation as

(I + K)ϕ = H 2P u inc + iηu inc where (I + K) = H(I + D ′ + iηS),
H is an integral operator with a strongly singular kernel which is injective and K is an integral operator with a weakly singular kernel. Uniqueness of a solution to the regularized integral equation follows from the uniqueness of a solution to (2.2). We conclude with Riesz theory and the fact that the equation (2.2) and the regularized one have the same solutions.

The far-field pattern can be computed via the integral representation formula

u ∞ = -F D P (u + u inc ),
where the far-field operator F D is defined by (see [2, equations (2.12) and (2.13)])

F D ϕ(x) = Γ 1 µ [I R 3 -x ⊗ x] e -iκs x•y 4π + 1 λ + 2µ [x ⊗ x] e -iκp x•y 4π ϕ(y) ds(y). (2.3)
Now we turn to the cavity problem: the total displacement field u + u inc is given by the superposition of the incident field u inc , which we assume to be an entire solution of the Navier equation, and the scattered field u, which solves the Navier equation in Ω c and satisfies the Neumann boundary condition (1.2) with g = -P u inc and the Kupradze radiation condition. A similar procedure as for the rigid body problem leads to the following integral representation of the solution

u(x) = Γ T [P y Φ(x, y)] u(y) + u inc (y) ds(y), x ∈ Ω c . ( 2
(D -I)(u + u inc ) = -2u inc and N (u + u inc ) = -2P u inc on Γ,
where D is the double layer boundary integral operator and N is the traction derivative of the double layer potential, defined by

D ϕ(x) = 2 Γ T [P y Φ(x, y)] ϕ(y) ds(y), N ϕ(x) = 2 Γ P x T [P y Φ(x, y)] ϕ(y) ds(y).
Setting ϕ = u + u inc , the field u given by (2.4) solves the Neumann boundary value problem (1.1)-( 1.3) if the density ϕ solves the following integral equation

N ϕ + iηD ϕ -iη ϕ = -2 P u inc + iηu inc on Γ. (2.5)
Here again, existence and uniqueness can be shown by the use of left or right regularization techniques. Indeed, the operator N :

H 1 2 (Γ) → H -1 2 (Γ)
is bounded and the operator D ′ -I, considered in these function spaces, is compact. By the identities SN = -I + D 2 and N S = -I + D ′ 2 , the single layer boundary integral operator for the Lamé system (ω = 0) denoted S 0 can be considered as a regularizer. Note that S 0 : Proposition 2.4] and S is a compact perturbation of S 0 . Therefore, the boundary integral equation operators became compact perturbations of the operators -I + D 2 or -I + D ′ 2 . Following [16, proof of Lemma 3.9] one can prove that these operators are Fredholm operators of index zero. Hence, we can apply Riesz theory.

H -1 2 (Γ) → H 1 2 (Γ) is invertible [1,
The far-field pattern can be computed via the integral representation formula

u ∞ = F N (u + u inc ),
where the far-field operator F N is defined by (see [2, equations (2.12) and (2.13)])

F N ϕ(x) = Γ 1 µ T P y [I R 3 -x ⊗ x] e -iκs x•y 4π + 1 λ + 2µ T P y [x ⊗ x] e -iκp x•y 4π ϕ(y) ds(y).
For both of the Dirichlet and the Neumann boundary value problems, the indirect approach is based on the following ansatz for the integral representation of the scattered elastic field :

u(x) = Γ T [P y Φ(x, y)]ϕ(y)ds(y) + iη Γ Φ(x, y)ϕ(y)ds(y), (2.6)
where η is still a non-zero real constant and ϕ is a density in H

N ϕ + iηD ′ ϕ -iηϕ = 2g on Γ. (2.8)
The far field is then given by u ∞ = F N ϕ + iηF D ϕ.

To apply the spectral method of Ganesh and Hawkins [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF], for solving one of the above CFIEs, we need to express the operators D, D ′ and N in terms of boundary integral operators with a weakly singular kernel and surface differential operators. Such representations can be obtained with the help of the tangential Günter derivative whose properties are given in the following lemma (see [START_REF] Kupradze | Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF]Equation (1.13) and Theorem 1.3, pages 282-284]) .

Lemma 2.1 Let Γ be a closed orientable surface in R 3 . The tangential Günter derivative denoted by M is defined for a vector function v ∈ C 1 (Γ, C 3 ) by

Mv = ∂ ∂n v -(div v)n + n × curl v.
For any vector functions v, ṽ in C 1 (Γ, C 3 ) there holds the Stokes formula

Γ (Mv) • ṽ ds = Γ v • (Mṽ) ds. (2.9)
The tangential gradient ∇ Γ and the surface divergence div Γ are defined for a scalar function u and a vector function v by the following equalities [START_REF] Nédélec | Acoustic and electromagnetic equations[END_REF]:

∇u = ∇ Γ u + ∂u ∂n n, div v = div Γ v + n • ∂v ∂n ,
and the Günter tangential derivative can then be rewritten as :

Mv = [∇ Γ v] -(div Γ v) • I R 3 n, (2.10)
where [∇ Γ v] is the matrix whose the j-th column is the tangential gradient of the j-th component of v.

In the following two lemmas we give the integral representation of the operators D, D ′ and N that we will consider in the remaining of the paper. We refer to [25, pages 47-49] for more details.

Lemma 2.2

The double layer boundary integral operator and the traction derivative of the single layer potential operator can be rewritten as :

(Dϕ)(x) = 2µ (SMϕ)(x) + 2 - Γ G(κ s , x -y)M y ϕ(y) ds(y) + Γ ∂ ∂n(y) G(κ s , x -y)ϕ(y)ds(y) + Γ ∇ y G(κ p , x -y) -G(κ s , x -y) n(y) • ϕ(y) ds(y) , (2.11)
and

(D ′ ϕ)(x) = 2µ (MSϕ)(x) + 2 -M x Γ G(κ s , x -y)ϕ(y) ds(y) + Γ ∂ ∂n(x) G(κ s , x -y)ϕ(y)ds(y) + n(x) Γ ϕ(y) • ∇ x G(κ p , x -y) -G(κ s , x -y) ds(y) .
(2.12)

Proof. We use integration by part and the properties (i) and (ii) of Lemma 2.1. The procedure to obtain these formulations is already known (we follow [25, page 49] in the case of the stationary Navier equations). We recall the proof since we use this development to prove the next lemma. First of all we rewrite the operator P as

P u = 2µMu + (λ + 2µ)(div u)n -µn × curl u. (2.13)
Then we apply the operator P y in the form (2.13) to the tensor Φ(x, y). It follows

T P y Φ(x, y) = 2µ T M y Φ(x, y) - T n(y) × curl y G(κ s , x -y)I R 3 + (λ + 2µ) µ T [n(y) ⊗ div y Φ(x, y)], div y Φ(x, y) = T ∇ y G(κ s , x -y) + 1 κ 2 s T ∇ y ∆ y G(κ s , x -y) -G(κ p , x -y) = κ 2 p κ 2 s T ∇ y G(κ p , x -y) , n(y) × curl y G(κ s , x -y)I R 3 = M y - ∂ ∂n(y) + n(y) • div y G(κ s , x -y)I R 3 .
In virtue of the property (i) in Lemma 2.1 we can write

T n(y) × curl y G(κ s , x -y)I R 3 = -M y - ∂ ∂n(y) G(κ s , x -y)I R 3 + T n(y) ⊗ T ∇ y G(κ s , x -y) .
Collecting the equalities we obtain

T P y Φ(x, y) = 2µ T M y Φ(κ, x -y) + ∂ ∂n(y) + M y G(κ s , x -y)I R 3 + ∇ y G(κ p , x -y) -G(κ s , x -y) ⊗ n(y).
To obtain (2.11) we use the equalities (2.9) and to obtain (2.12) we take the adjoint form of (2.11).

Lemma 2.3

The traction derivative of the double layer potential operator can be rewritten as follows:

N = 2µD ′ M + 2µM(D -2µSM) + V with V ϕ(x) = ρω 2 Γ n(x)G(κ p , x -y) n(y) • ϕ(y) ds(y) + ρω 2 n(x) × Γ G(κ s , x -y)(ϕ(y) × n(y))ds(y) -µ curl Γ Γ G(κ s , x -y) curl Γ ϕ(y)ds(y), x ∈ Γ, (2.14)
where we have

curl Γ = -n × ∇ Γ and curl Γ = div Γ ( • × n) .
Proof. First, we write D = 2µ SM + (D -2µ SM). Applying the traction trace to the kernel of 2µ SM we obtain 2µ D ′ M. Then, we apply the operator P x , in the form (2.13), to the kernel of (D -2µ SM). Following the calculation of the previous proof, one can see that we have

(D -2µ SM)ϕ(x) = - Γ T [n(y) × curl y {G(κ s , x -y)I 3 }]ϕ(y)ds(y) + Γ ∇ y G(κ p , x -y) n(y) • ϕ(y) ds(y).
We observe that

T [n(y) × curl y {G(κ s , x -y)I 3 }]ϕ(y) = T [curl y {G(κ s , x -y)I 3 }](ϕ(y) × n(y)) = -[curl y {G(κ s , x -y)I 3 }] (ϕ(y) × n(y)) = [curl x {G(κ s , x -y)I 3 }](ϕ(y) × n(y)) = curl x {G(κ s , x -y)(ϕ(y) × n(y))} .
The composition of the operator (λ + 2µ)n(x) div x with this term vanishes since one has div curl = 0 and the composition with µn(x) × curl x yields

µn(x) × (-∆ x + ∇ x div x ) {G(κ s , x -y)(ϕ(y) × n(y))} = µκ 2 s n(x) × {G(κ s , x -y)(ϕ(y) × n(y))} -µ curl x Γ div x {G(κ s , x -y)(ϕ(y) × n(y))} = ρω 2 n(x) × {G(κ s , x -y)(ϕ(y) × n(y))} -µ curl x Γ {∇ x G(κ s , x -y) • (ϕ(y) × n(y))} = ρω 2 n(x) × {G(κ s , x -y)(ϕ(y) × n(y))} -µ curl x Γ {-∇ y G(κ s , x -y) • (ϕ(y) × n(y))} = ρω 2 n(x) × {G(κ s , x -y)(ϕ(y) × n(y))} -µ curl x Γ {-∇ y Γ G(κ s , x -y) • (ϕ(y) × n(y))} , (2.15)
since

curl x curl x = -∆ x + ∇ x div x and ∇ Γ • = n × (∇ • ×n).
In view of the following integral part formula for a scalar density ϕ 1 and a tangential vector density ϕ 2

- Γ (∇ Γ ϕ 1 ) • ϕ 2 ds = + Γ ϕ 1 div Γ ϕ 2 ds,
and the identity div Γ (ϕ 2 × n) = curl Γ ϕ 2 , we obtain the last two terms of the right hand side in (2.14).

We observe that ∇ y G(κ p , x-y) n(y)•ϕ(y) = -∇ x G(κ p , x-y) n(y)•ϕ(y) . The composition of the operator µn(x) × curl x with this term vanishes since one has curl ∇ = 0 and the composition with (λ + 2µ)n(x) div x yields

-(λ + 2µ)n(x) div x ∇ x G(κ p , x -y) n(y) • ϕ(y) = -(λ + 2µ)n(x)∆ x G(κ p , x -y) n(y) • ϕ(y) = ρω 2 n(x) G(κ p , x -y) n(y) • ϕ(y) .
We conclude by collecting all these formulas.

Spherical parametrization of the integral equations

The spectrally accurate algorithm is based on a spherical reformulation of the integral equations. We denote by S 2 the unit sphere of R 3 . Let q : S 2 → Γ be a parametrization of class C 1 at least. We denote by τ q the transformation ("pull-back") that maps a function u defined on Γ onto the function τ q (u) = u • q defined on S 2 . Let J q be the determinant of the Jacobian of the change of variable x ∈ S 2 → q(x) ∈ Γ. We denote by θ, φ the spherical coordinates of any point x ∈ S 2 , i.e.

x = ψ(θ, φ) =   sin θ cos φ sin θ sin φ cos θ   , (θ, φ) ∈ ]0; π[×[0; 2π[ ∪{(0, 0); (0, π)}.
The tangent and the cotangent planes at any point x = ψ(θ, φ) ∈ S 2 are generated by the unit vectors :

e θ = ∂ψ ∂θ (θ, φ) =   cos θ cos φ cos θ sin φ -sin θ   and e φ = 1 sin θ ∂ψ ∂φ (θ, φ) =   -sin φ cos φ 0   .
The triplet (x, e θ , e φ ) forms a direct orthonormal system. The determinant of the Jacobian of the change of variable is J ψ (θ, φ) = sin θ.

The matrix [D S 2 q(x)] = T [∇ S 2 q(x)] maps the tangent plane T x to S 2 at the point x onto the tangent plane T q(x) to Γ at the point q(x). The latter is generated by the vectors

t 1 (x) = ∂q • ψ ∂θ • ψ -1 and t 2 (x) = 1 sin θ ∂q • ψ ∂φ • ψ -1 . (3.1)
We can write [D S 2 q(x)] = t 1 ⊗ e θ + t 2 ⊗ e φ , where for two vectors a, b ∈ R 3 we set a ⊗ b = a T b. The functions J q and τ q (n) can be computed via the formulas : J q = t 1 × t 2 and τ q (n) = (t 1 × t 2 )/J q . The parametrization q : S 2 → Γ being a diffeomorphism, we set [D S 2 q(x)] -1 = [D Γ q -1 ] • q(x). The transposed matrix T [D S 2 q(x)] -1 maps the cotangent plane T * x to S 2 at the point x onto the cotangent plane T * q(x) to Γ at the point q(x). The latter is generated by the vectors

t 1 (x) = t 2 (x) × n(q(x)) J q (x) and t 2 (x) = n(q(x)) × t 1 (x) J q (x) . (3.2)
We can write

[D S 2 q(x)] -1 = e θ ⊗ t 1 + e φ ⊗ t 2 .
In addition to the formulas [23, Eq. (3.

2)] we will use, for non-tangential vector functions v, w ∈ C 1 (Γ, C 3 ), the following identities

τ q (div Γ v) = [D S 2 q] -1 : [∇ S 2 τ q (v)], (3.3)
τ q (curl Γ w) = - 1 J q T [D S 2 q] : curl S 2 τ q (w) , (3.4

Spherical reformulation

The key of the spectral method is to transport the boundary integral equation on the unit sphere. To do so, following [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF][START_REF] Ganesh | A spectrally accurate algorithm for electromagnetic scattering in three dimensions[END_REF][START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF], we use a change of variable in the integrand and we consider instead the following integral operators:

S ϕ s = J q τ q Sτ -1 q (ϕ s ) = 2 J q S 2 Φ(q(•), q(ŷ)) ϕ s (ŷ) J q (ŷ) ds(ŷ), D ϕ s = J q τ q Dτ -1 q (ϕ s ) = 2µ SM q ϕ s + 2 J q - S 2
G κ s , q(•)q(ŷ) M q(ŷ) ϕ s (ŷ) J q (ŷ)ds(ŷ)

+ S 2
∂ ∂n(q(ŷ)) G(κ s , q(•)q(ŷ))ϕ s (ŷ) J q (ŷ)ds(ŷ)

+ S 2 ∇ q(ŷ) G(κ p , q(•) -q(ŷ)) -G(κ s , q(•) -q(ŷ)) n(q(ŷ)) • ϕ s (ŷ) J q (ŷ)ds(ŷ) , D ′ ϕ s = J q τ q D ′ τ -1 q (ϕ s ) = 2µM q S ϕ s + 2 J q -M q S 2 G κ s , q(•) -q(ŷ) ϕ s (ŷ) J q (ŷ)ds(ŷ) + S 2 ∂ ∂n(q(•)) G(κ s , q(•) -q(ŷ))ϕ s (ŷ) J q (ŷ)ds(ŷ) + n(q(•)) S 2 ϕ s (q(ŷ)) • ∇ q(•) G(κ p , q(•) -q(ŷ)) -G(κ s , q(•) -q(ŷ)) J q (ŷ)ds(ŷ) , and N ϕ s = J q τ q N τ -1 q ϕ s ) = 2µD ′ M q ϕ s + 2µM q (D -2µSM q )ϕ s + Vϕ s with Vϕ s = J q τ q N 2µ τ -1 q ϕ s ) = 2 J q ρω 2 n(q(•)) S 2 G(κ p , q(•) -q(ŷ))(n(q(ŷ)) • ϕ s (ŷ)) J q (ŷ)ds(ŷ) + n(q(•)) × S 2 G κ s , q(•) -q(ŷ) ϕ s (ŷ) × n(q(ŷ)) J q (ŷ)ds(ŷ) + µ[D S 2 q] curl S 2 S 2 G(κ s , q(•) -q(ŷ)) T [D S 2 q(ŷ)] : curl S 2 ϕ s (ŷ)] ds(ŷ), where M q(ŷ) ϕ s (ŷ) = T [D S 2 q(ŷ)] -1 [∇ S 2 ϕ s (ŷ)] -[D S 2 q(ŷ)] -1 : [∇ S 2 ϕ s (ŷ)] I R 3 n(q(ŷ)). (3.5)
Splitting of the kernels Following details in [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF][START_REF] Ganesh | A spectrally accurate algorithm for electromagnetic scattering in three dimensions[END_REF][START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF], to implement the integral operators S, D, D ′ and N we split their kernels into a smooth and a weakly singular part. We introduce the functions :

G 1 (q; κ, x, ŷ) = 1 2π cos(κ|q(x) -q(ŷ)|), G 2 (q; κ, x, ŷ) = 1 2π    sin(κ|q(x) -q(ŷ)|) |q(x) -q(ŷ)| if x = ŷ, κ if x = ŷ. and R(q; x, ŷ) = |x -ŷ| |q(x) -q(ŷ)| .
Then S can be rewritten as:

S ϕ s (x) = S 2 R(q; x, ŷ) |x -ŷ| W 1 (q; x, ŷ)ϕ s (ŷ)ds(ŷ) + i S 2 W 2 (q; x, ŷ)ϕ s (ŷ)ds(ŷ)
where W 1 (q; x, ŷ) and W 2 (q; x, ŷ) are 3 × 3 matrix-valued functions given by

W 1 (q; x, ŷ) = J q (x)J q (ŷ) µ G 1 (q; κ s , x, ŷ) - 1 κ 2 s G 1 (q; κ s , x, ŷ) -G 1 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 2 I R 3 - J q (x)J q (ŷ) µκ 2 s κ s G 2 (q; κ s , x, ŷ) -κ p G 2 (q; κ p , x, ŷ) I R 3 - J q (x)J q (ŷ) µκ 2 s κ 2 s G 1 (q; κ s , x, ŷ) -κ 2 p G 1 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 2 A(q, x, ŷ) +3 J q (x)J q (ŷ) µκ 2 s G 1 (q; κ s , x, ŷ) -G 1 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 4 A(q, x, ŷ) +3 J q (x)J q (ŷ) µκ 2 s κ s G 2 (q; κ s , x, ŷ) -κ p G 2 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 2 A(q, x, ŷ), W 2 (q; x, ŷ) = J q (x)J q (ŷ) µ G 1 (q; κ s , x, ŷ) - 1 κ 2 s G 2 (q; κ s , x, ŷ) -G 2 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 2 I R 3 + J q (x)J q (ŷ) µκ 2 s κ s G 1 (q; κ s , x, ŷ) -κ p G 1 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 2 I R 3 -3 J q (x)J q (ŷ) µκ 2 s κ s G 1 (q; κ s , x, ŷ) -κ p G 1 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 4 A(q; x, ŷ) - J q (x)J q (ŷ) µκ 2 s κ 2 s G 2 (q; κ s , x, ŷ) -κ 2 p G 2 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 2 A(q; x, ŷ) +3 J q (x)J q (ŷ) µκ 2 s G 2 (q; κ s , x, ŷ) -G 2 (q; κ p , x, ŷ) |q(x) -q(ŷ)| 4 A(q; x, ŷ), with A(q; x, ŷ) = q(x) -q(ŷ) ⊗ q(x) -q(ŷ) .
The operator D can be rewritten as:

D ϕ s (x) = S 2 R(q; x, ŷ) |x -ŷ| 2µW 1 (q; x, ŷ) -G 1 (q; κ s , x, ŷ)J q (x)J q (ŷ)I R 3 M q(ŷ) ϕ s (ŷ)ds(ŷ) + i S 2 2µW 2 (q; x, ŷ) -G 2 (q; κ s , x, ŷ)J q (x)J q (ŷ)I R 3 M q(ŷ) ϕ s (ŷ)ds(ŷ) + S 2 R(q; x, ŷ) |x -ŷ| D 1 (q; x, ŷ) + K 1 (q; x, ŷ) I R 3 ϕ s (ŷ)ds(ŷ) + i S 2 D 2 (q; x, ŷ) + K 2 (q; x, ŷ) I R 3 ϕ s (ŷ)ds(ŷ)
where K 1 (q; x, ŷ) and K 2 (q; x, ŷ) are scalar-valued functions given by K 1 (q; x, ŷ) = J q (x)J q (ŷ) G 1 (q; κ s , x, ŷ) |q(x)q(ŷ)| 2 + κ s G 2 (q; κ s , x, ŷ) n(q(ŷ)) • q(x)q(ŷ)

K 2 (q; x, ŷ) = J q (x)J q (ŷ) G 2 (q; κ s , x, ŷ) -κ s G 1 (q; κ s , x, ŷ) |q(x) -q(ŷ)| 2 n(q(ŷ)) • q(x) -q(ŷ)
and D 1 (q; x, ŷ) and D 2 (q; x, ŷ) are 3 × 3 matrix-valued functions given by

D 1 (q; x, ŷ) = J q (x)J q (ŷ) G 1 (q; κ p , x, ŷ) -G 1 (q; κ s , x, ŷ) |q(x) -q(ŷ)| 2 q(x) -q(ŷ) ⊗ n(q(ŷ)) + J q (x)J q (ŷ) κ p G 2 (q; κ p , x, ŷ) -κ s G 2 (q; κ s , x, ŷ) q(x) -q(ŷ) ⊗ n(q(ŷ)), D 2 (q; x, ŷ) = J q (x)J q (ŷ) G 2 (q; κ p , x, ŷ) -G 2 (q; κ s , x, ŷ) |q(x) -q(ŷ)| 2 q(x) -q(ŷ) ⊗ n(q(ŷ)) -J q (x)J q (ŷ) κ p G 1 (q; κ p , x, ŷ) -κ s G 1 (q; κ s , x, ŷ) |q(x) -q(ŷ)| 2 q(x) -q(ŷ) ⊗ n(q(ŷ)).
The operator D ′ can be rewritten as:

D ′ ϕ s (x) = M q(x) S 2 R(q; x, ŷ) |x -ŷ| 2µW 1 (q; x, ŷ) -G 1 (q; κ s , x, ŷ)J q (x)J q (ŷ)I R 3 ϕ s (ŷ)ds(ŷ) + i M q(x) S 2 2µW 2 (q; x, ŷ) -G 2 (q; κ s , x, ŷ)J q (x)J q (ŷ)I R 3 ϕ s (ŷ)ds(ŷ) + S 2 R(q; x, ŷ) |x -ŷ| D ′ 1 (q; x, ŷ) + K ′ 1 (q; x, ŷ) I R 3 ϕ s (ŷ)ds(ŷ) + i S 2 D ′ 2 (q; x, ŷ) + K ′ 2 (q; x, ŷ) I R 3 ϕ s (ŷ)ds(ŷ)
where K ′ 1 (q; x, ŷ) and K ′ 2 (q; x, ŷ) are scalar-valued functions defined by

K ′ 1 (q; x, ŷ) = K 1 (q; ŷ, x), K ′ 2 (q; x, ŷ) = K 2 (q; ŷ, x).
and D ′ 1 (q; x, ŷ) and D ′ 2 (q; x, ŷ) are 3 × 3 matrix-valued functions defined by

D ′ 1 (q; x, ŷ) = T D 1 (q; ŷ, x), D ′ 2 (q; x, ŷ) = T D 2 (q; ŷ, x).
The operator V can be rewritten as:

V ϕ s (x) = S 2 R(q; x, ŷ) |x -ŷ| V 1 (q; x, ŷ)ϕ s (ŷ)ds(ŷ) + i S 2 V 2 (q; x, ŷ)ϕ s (ŷ)ds(ŷ) + µ [D S 2 q(x)] curl S 2 S 2 R(q; x, ŷ) |x -ŷ| G 1 (q; κ s , x, ŷ) T [D S 2 q(ŷ)] : [curl S 2 ϕ s (ŷ)] ds(ŷ) + iµ [D S 2 q(x)] curl S 2 S 2 G 2 (q; κ s , x, ŷ) T [D S 2 q(ŷ)] : [curl S 2 ϕ s (ŷ)] ds(ŷ),
where V 1 (q; x, ŷ) and V 2 (q; x, ŷ) are 3 × 3 matrix-valued functions defined by V 1 (q; x, ŷ) = ρω 2 G 1 (q; κ p , x, ŷ)J q (x)J q (ŷ) n(q(x)) ⊗ n(q(ŷ))

+ G 1 (q; κ s , x, ŷ) J q (x)J q (ŷ) n(q(x)) • n(q(ŷ)) I R 3n(q(ŷ)) ⊗ n(q(x)) , V 2 (q; x, ŷ) = ρω 2 G 2 (q; κ p , x, ŷ)J q (x)J q (ŷ) n(q(x)) ⊗ n(q(ŷ))

+ G 2 (q; κ s , x, ŷ) J q (x)J q (ŷ) n(q(x)) • n(q(ŷ)) I R 3n(q(ŷ)) ⊗ n(q(x)) .

Finally, setting ϕ s = P (u + u inc ) • q and 2f s = J q (P u inc + iηu inc ) • q, the parametrized form of equation (2.2) is

J q I H 1 2 (S 2 ) + D ′ + iηS ϕ s = f s , on S 2 , (3.6)
and setting ϕ s = (u + u inc ) • q and g s = -2J q (P u inc + iηu inc ) • q the parametrized form of equation (2.5) is

N + iηD -iη J q I H 1 2 (S 2 ) ϕ s = g s , on S 2 , (3.7)
The parametrized forms of the equations (2.7) and (2.8) are, respectively, the adjoint forms of the equations (3.6) and (3.7).

The high-order spectral algorithm

The spectral algorithm is based on that in [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF]. In a first step, we deal with the singularities of the weakly singular parts of the integrands by introducing a change of coordinate system on S 2 . It will yield transformed operators with kernels which are singular only at one point on the sphere, namely the north pole. In a second step, we interpolate the integrands and project the integral equations on vector spherical function basis.

Treatment of the singularities For any x ∈ S 2 , let T x be the orthogonal transformation defined in [17, Subsection 3.1.1] that satisfies T x x = T (0, 0, 1) = η. For any ŷ ∈ S 2 , we set ẑ = T x ŷ. We also introduce an induced linear tranformation T x defined by T xu(ẑ) = u(T -1

x ẑ) = u(ŷ) and we still denote by T x its bivariate analogue T xv(ẑ 1 , ẑ2 ) = v(T -1

x ẑ1 , T -1 x ẑ2 ). The boundary integral operator S can be rewritten in the form:

S ϕ s (x) = S 2 T x R(q; η, ẑ) |η -ẑ| T x W 1 (q; η, ẑ)T x ϕ s (ẑ)ds(ẑ) + i S 2 T x W 2 (q; η, ẑ)T x ϕ s (ẑ)ds(ẑ), (4.1)
and it can be shown that

(θ ′ , φ ′ ) → T x R(q; η, ẑ(θ ′ , φ ′ ))T x W 1 (q; η, ẑ(θ ′ , φ ′ )) is smooth. An important point is that the singularity 1 |η -ẑ(θ ′ , φ ′ )| = 1 2 sin θ ′ 2
is cancelled out by the surface element ds(ẑ) = sin θ ′ dθ ′ dφ ′ . We proceed in the same way for the integrand of the operators D, D ′ and V.

Fully discrete approximations Let O n and L n be respectively the operators defined in [20, Equations (3.10) and (3.12)]. For practical purposes, we choose n ′ = 2n + 1 and n ≥ 5. (This satisfies the theoretical constraints required on n and n ′ for convergence analysis based on that in [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF][START_REF] Ganesh | A spectrally accurate algorithm for electromagnetic scattering in three dimensions[END_REF][START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF].) The codomain H n of the fully discrete operator O n is the 3(n + 1) 2 -2 finite dimensional space spanned by orthonormal vector spherical harmonics of degree at most n. Following details in [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF], the spectral algorithm for the rigid body problem is: find

ϕ n ∈ H n such that O n ϕ n + O n D ′ n ′ ϕ n + iηO n S n ′ ϕ n = 2O n f s , (4.2)
where, as in [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF], S n ′ and D ′ n ′ are defined as

S n ′ ϕ s (x) = S 2 1 |η -ẑ| L n ′ {T x R(q; η, •)T x W 1 (q; η, •)T x ϕ s (•)}(ẑ) ds(ẑ) + i S 2 L n ′ {T x W 2 (q; η, •)T x ϕ s (•)}(ẑ) ds(ẑ), and
D ′ n ′ ϕ s (x) = S 2 1 |η -ẑ| L n ′ T x R(q; η, •)T x D ′ 1 (q; η, •) + K ′ 1 (q; η, •) I R 3 T x ϕ s (•) (ẑ) ds(ẑ) + i S 2 L n ′ T x D ′ 2 (q; η, •) + K ′ 2 (q; η, •) I R 3 T x ϕ s (•) (ẑ)ds(ẑ) + M q(x) S 2 1 |η -ẑ| L n ′ T x R(q; η, •)T x 2µW 1 (q; η, •) -G 1 (q; κ s , η, •)J q (•)I R 3 T x ϕ s (•) (ẑ) ds(ẑ) + i M q(x) S 2 L n ′ T x 2µW 2 (q; η, •) -G 2 (q; κ s , η, •)J q (•)I R 3 T x ϕ s (•) (ẑ)ds(ẑ).
The spectral algorithm for the cavity problem is: find

ϕ n ∈ H n such that O n N n ′ ϕ n + iηO n D n ′ ϕ n -iηO n ϕ n = 2O n g s , (4.3)
where

D n ′ ϕ s (x) = S 2 1 |η -ẑ| L n ′ T x R(q; η, •)T x D 1 (q; η, •) + K 1 (q; η, •) I R 3 T x ϕ s (•) (ẑ) ds(ẑ) + i S 2 L n ′ T x D 2 (q; η, •) + K 2 (q; η, •) I R 3 T x ϕ s (•) (ẑ) ds(ẑ) + S 2 1 |η -ẑ| L n ′ T x R(q; η, •)T x 2µW 1 (q; η, •) -G 1 (q; κ s , η, •)J q (•)I R 3 T x (M q ϕ s )(•) (ẑ) ds(ẑ) + i S 2 L n ′ T x 2µW 2 (q; η, •) -G 2 (q; κ s , η, •)J q (•)I R 3 T x (M q ϕ s )(•) (ẑ) ds(ẑ),
and

N n ′ = 2µ(D ′ M) n ′ + 2µM(D -2µSM) n ′ + V n ′
denotes the corresponding approximation of the operator N with

V n ′ ϕ s (x) = S 2 1 |η -ẑ| L n ′ T x R(q; η, •)T x V 1 (q; η, •)T x ϕ s (•) (ẑ) ds(ẑ) + i S 2 L n ′ T x V 2 (q; η, •)T x ϕ s (•) (ẑ) ds(ẑ) + µ [D S 2 q(x)] curl S 2 S 2 1 |η -ẑ| L n ′ T x R(q; η, •)T x G 1 (q; κ s , η, •)T x T [D S 2 q] : [curl S 2 ϕ s ] (•) (ẑ) ds(ẑ) + iµ [D S 2 q(x)] curl S 2 S 2 L n ′ T x G 2 (q; κ s , η, •)T x T [D S 2 q] : [curl S 2 ϕ s ] (•) (ẑ) ds(ẑ).
As in [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF], the above discrete operators can be written using summation form using the eigenfunction properties of the spherical harmonics [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]. For details of convergence analysis of the method developped originally for electromagnetic scattering, we refer to [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF][START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF] where such bounds are established.

Remark 4.1 Contrary to the acoustic and the electromagnetic cases, the algorithm additionnally requires the evaluation of the tangential Günter derivative M q (3.5) and the tangential vector curl (3.4) of both the approximate densities ϕ n in the integrands (i.e. with respect to the variable ŷ) and some densities resulting from surface integrals (i.e. with respect to the variable x).

(ii) A numerical evaluation of the tangential derivatives with respect to x is obtained through the application of the discrete projection operator O n and the integration by parts (2.9) and

Γ curl Γ u • ṽ ds = Γ u curl Γ ṽ ds,
put in an other way

S 2 M q(x) v(x) • ṽ(x) J q (x) ds(x) = S 2 v(x) • M q(x) ṽ(x) J q (x) ds(x) S 2 [D S 2 q(x)] curl S 2 u(x) • ṽ(x) ds(x) = - S 2 u(x) T [D S 2 q(x)] : [curl S 2 ṽ(x)] ds(x) .
Thus, we are lead to the evaluation of the tangential derivatives of the spherical vector harmonics. These integration by parts formulas justify why we multiply the boundary integral equations by J q .

Numerical experiments

In the remaining of the paper we present numerical experiments for low and medium frequencies to highlight the fast convergence of our algorithm through the evaluation of the far-field pattern u ∞ and the boundary traces of the scattered field. We use the various smooth and non-smooth, convex and non convex obstacles whose the parametric representations and visualizations are given in [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF]. The surfaces are characterized by their diameter denoted size obj. The convex shapes are the sphere, denoted by sphere(size obj) and the ellipsoid with radius a, b and c denoted by ellipsoid(a, b, c). The non convex shapes are the bean denoted by bean(size obj) and the peanut denoted by pea(size obj, α). As a first test, using the CFIEs based on the indirect approach, we compute the far field, denoted by u ∞ ps , created by an off center point source located inside the elastic obstacle :

u inc (x) = -[Φ(x, s)]p, s ∈ Ω and p ∈ S 2 .
In this case the total exterior wave has to vanish so that the far-field pattern of the scattered wave u s is the opposite of the far field pattern of the incident wave. The following far field representation is obtain by applying the kernel of the far field operator F D to p (see (2.3)).

u ∞ exact (x) = 1 µ e -iκs x•s 4π (x × p) × x + 1 λ + 2µ e -iκp x•s 4π (x • p) x.
We choose s = T (0, 0.05, 0.0866) and p = T (1, 0, 0). In the tabulated results we indicate the uniform-norm error (by taking the maximum of errors obtained over 1300 observed directions, i.e n ∞ = 25) :

||[u ∞ ps ] n -u ∞ exact || ∞ = max x∈S 2 [u ∞ ps ] n -u ∞ exact .
As a second test, using CFIEs based on the direct approach, we compute the far field, denoted by u ∞ pw , created by the scattering of an incident plane elastic wave [2, Section 3] or [24, section 3]:

u inc (x) = 1 µ e iκsx•d (d × p) × d + 1 λ + 2µ e iκpx•d (d • p)d, where d, p ∈ S 2 . (5.1)
When the polarization p is orthogonal to the propagation vector d, the incident plane waves oscillate in a direction orthogonal to the direction of propagation. They are called shearing waves. When p = d, the incident plane waves oscillate along the direction of propagation. They are called pressure waves. In the tabulated results (except for the sphere) we indicate the real part and the imaginary part of the polarization component of the far field evaluated at the incident direction : [u ∞ pw (d)] n • p. In the case of the sphere we know the analytical representation of u ∞ pw (d), thus we compute the errors ||[u ∞ pw ] nu ∞ exact || ∞ by taking the maximum errors over 1300 directions as for the point source radiation. The representation of u ∞ pw (d) is given by a series expansion of the vector spherical harmonics (see the Appendix for more details). Here, we consider an incident shearing plave wave defined with d = T (0, 0, 1) and p = T (1, 0, 0). In all simulations, we set η = 1 in the CFIE formulations and the mechanical parameters are chosen such that the wavenumbers satisfy κ s = 2κ p . The size of the scatterers is indicated between brackets and is expressed in terms of the S-wavelenght λ s = 2π/κ s . The experiments are realised with Matlab using a direct solver for the matrix inversion and parallel for-loops for the matrix setup. The dimension of the largest system we have to solve here is about 25 000 (n = 90). The CPU time is mainly used for the matrix setup. In the low frequency region (with S-wavelengths greater than 1/10 times the size of obstacles), the simulation of an accurate approximate solution takes from few seconds to few minutes of CPU time (on a laptop with 8 GB RAM, 2.9 GHz Intel core i7 processor with 4 threads in our case) and in about few hours of CPU time for high frequency scattering (with S-wavelengths close to 1/25 times the size of obstacles) using a 2.4 GHz Intel Xeon processor with 8 threads and 48 Go RAM. The results are presented in tables 1 to 3 for the Dirichlet problem and in tables 5 to 7 for the Neumann problem. As expected, we obtain similar numerical convergence rates than Ganesh & al in [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF][START_REF] Ganesh | A spectrally accurate algorithm for electromagnetic scattering in three dimensions[END_REF][START_REF]An efficient surface integral equation method for the time-harmonic Maxwell equations[END_REF][START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF][START_REF]A high-order tangential basis algorithm for electromagnetic scattering by curved surfaces[END_REF]. We observe that the presence of angularity or concavity greatly increase the needed number of degree of freedom.

As a third and last test, we compare our algorithm (FSM) with the FMM algorithm [START_REF] Chaillat | A fast multipole method formulation for 3d elastodynamics in the frequency domain[END_REF][START_REF] Chaillat | A multi-level fast multipole bem for 3-d elastodynamics in the frequency domain[END_REF][START_REF] Chaillat | Fast Multipole Method for 3-D elastodynamic boundary integral equations[END_REF] in the low frequency region. We consider the scattering problem of a shearing plane wave, with the above-detailed properties, by a rigid elastic sphere. In table 4 we give the relative L 2 -norms, denoted by ǫ 2 , of the traction trace and of the far-field of the scattered field in accordance with the number of degree of freedom. The data n λs indicates the density of discretization points per S-wavelength used for the FMM-BEM method. The FMM runned with an 3.49 GHz Intel Xeon processor and 32GB RAM. The results demonstrate that the spectral algorithm still is very competitive with the fast multipole method in linear elasticity. Indeed, the first one requires significantly fewer unknowns than the second one to achieve the same precision. This induces a total CPU time allocated to the matrix inversion (or iteration procedure) greatly lower with the FSM than with the FMM. Even if the iteration procedure of the FMM can be speed up using preconditionners ( work in progress based on [START_REF] Darbas | Analytic preconditioners for the iterative solution of elastic scattering problems[END_REF]), one can improve the full matrix setup with the FSM by performing parallel computations and using other fast algorithm as the Fast Fourier Transform or the Fast Legendre Transform.

= 24π n ||[u ∞ ps ] n -u ∞ exact || ∞ ||[u ∞ pw ] n -u ∞ exact || ∞ n ||[u ∞ ps ] n -u ∞ exact || ∞ ||[u ∞ pw ] n -u ∞ exact || ∞ 30
= 8π n ||[u ∞ ps ] n -u ∞ exact || ∞ Re[u ∞ pw (d)] n • p Im[u ∞ pw (d)] n • p 25
) n ||[u ∞ ps ] n -u ∞ exact || ∞ ||[u ∞ ps ] n -u ∞ exact || ∞ ||[u ∞ ps ] n -u ∞
= 24π n ||[u ∞ ps ] n -u ∞ exact || ∞ ||[u ∞ pw ] n -u ∞ exact || ∞ n ||[u ∞ ps ] n -u ∞ exact || ∞ ||[u ∞ pw ] n -u ∞ exact || ∞ 30

Conclusion and perspectives

In this paper we have proposed a new high order spectral algorithm, based on previous investigations in acoustics and electromagnetics [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF][START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF], for solving scattering problems of time-harmonic elastic waves in three dimensions. The approach has been successfully validated for a variety of geometries. A numerical comparison attests that the proposed fast spectral algorithm is more efficient than already existing fast BEM-based algorithm in linear elasticity. However, we point out that our algorithm only applies to scatterers with a smooth and simply connected closed boundary and requires the knowledge of an explicit spherical parametrisation of the boundary. These constraints limits the potential industrial applications of the fast spectral algorithm. Unbounded domains are often encountered in scientific or engineering applications. Typical examples in computational elastodynamics include earthquake simulation, remote sensing, and scattering from cracks or inhomogeneities. To solve an equation in an infinite domain numerically, it is usual to limit the computation to a finite domain Ω with a smooth and simply connected closed boundary Γ. The scattering problem can be reduced to the solution of the elastodynamic system in the interior of Ω coupled with an integral equation posed on the artificial boundary Γ which could be solved with such a spectral algorithm.

In a forthcoming paper, this algorithm will be used to develop and analyze novel Newton-type iterative methods, typically based on parametrizations of the unknown boundary [START_REF] Charalambopoulos | On the Fréchet differentiability of boundary integral operators in the inverse elastic scattering problem[END_REF][START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF][START_REF] Ivanyshyn | Nonlinear integral equation methods for the reconstruction of an acoustically sound-soft obstacle[END_REF], for solving inverse obstacle scattering problems in linear elasticity. An other extension of this work is to consider the scattering problem by multiple obstacles for various boundary conditions, with the objective to incorporate it in an inverse scattering algorithm (see [START_REF] Carpio | Determining planar multiple sound-soft obstacles from scattered acoustic fields[END_REF] for example).

. 1 )

 1 Taking the Neumann and Dirichlet traces of the above identity we obtain [15, formulas (2.8)-(2.10)]

. 4 )

 4 Taking the Neumann and Dirichlet traces of the above identity we obtain [15, formulas (2.8)-(2.10)]

1 2 (

 2 Γ). Taking the Dirichlet trace to (2.6) together with the boundary condition (1.2) we obtain the following CFIE for the Dirichlet boundary value problem ϕ + Dϕ + iηSϕ = 2f on Γ. (2.7) Taking the Neumann trace to (2.6) together with the boundary condition (1.3) we obtain the following CFIE for the Neumann boundary value problem

  ) where for two (3×3) matrices A and B whose columns are denoted by (a 1 , a 2 , a 3 ) and (b 1 , b 2 , b 3 ), respectively, we set A : B = a 1 • b 1 + a 2 • b 2 + a 3 • b 3 . These equalities are an extension of the formulas stated in [23, Equation (3.2) and Appendix A] by using the well known indentities div Γ v = T race([∇ Γ v]) and curl Γ w = -T race([curl Γ w]).

Table 1 :

 1 Scattering by a sphere (Dirichlet Problem). sphere(8.0λ s ) , κ s = 8π sphere(24.0λ s ) , κ s

Table 2 :

 2 Scattering by an ellipsoid (Dirichlet Problem). ellipsoid(8.0λ s , 6.0λ s , 4.0λ s ) , κ s

		1.9135E -14	6.1749E -03	80	2.2754E -12	5.2605E -01
	35	2.0324E -14	3.2329E -06	85	3.2748E -12	2.1300E -03
	40	2.0139E -14	1.4363E -08	90	3.0681E -12	3.8588E -06

Table 3 :

 3 Scattering by various boundaries (Dirichlet Problem).

		4.2206E -03	-1.562 106 455	10.515 343 810
	30	2.3785E -05	-1.564 562 097	10.516 585 108
	35	7.7785E -09	-1.564 570 704	10.516 577 444
	40	6.0542E -12	-1.564 570 705	10.516 577 443

Table 4 :

 4 Comparison of the present spectral algorithm (FSM) with the Fast Multipole method (FMM)[START_REF] Chaillat | A fast multipole method formulation for 3d elastodynamics in the frequency domain[END_REF]. of an incident S-plane wave by a rigid sphere for various frequencies.

	exact || ∞

Table 5 :

 5 Scattering by a sphere (Neumann Problem). sphere(8.0λ s ) , κ s = 8π sphere(24.0λ s ) , κ s

Table 6 :

 6 Scattering by an ellipsoid (Neumann Problem). ellipsoid(8.0λ s ,6λ s ,4.0λ s ) ,κ s = 8π n ||[u ∞ ps ] nu ∞ exact || ∞

		2.0938E -13	1.3201E -02	80	4.1592E -12	8.0675E -01
	35	4.6209E -14	6.5179E -06	85	6.5344E -12	2.8372E -03
	40	2.1012E -13	1.0740E -08	90	5.1414E -12	3.1626E -06
				Re[u ∞ pw (d)] n • p	Im[u ∞ pw (d)] n • p
	25		5.6889E -03	1.401 764 629	10.096 850 244
	30		2.5794E -04	1.294 233 563	10.198 853 416
	35		5.7971E -08	1.293 643 766	10.192 397 631
	40		1.4182E -10	1.293 643 351	10.192 399 357

Table 7 :

 7 Scattering by various boundaries (Neumann Problem). κ s = 24π κ s = 16π ellipsoid(24.0λ s ,18.0λ s ,12.0λ s ) pea(16.0λ s ,0.25) bean(16.0λ s )

	n	||[u ∞ ps ] n -u ∞ exact || ∞	||[u ∞ ps ] n -u ∞ exact || ∞	||[u ∞ ps ] n -u ∞ exact || ∞
	80	4.8694E -05	3.6309E -01	3.1169E -03
	85	5.8227E -06	9.0615E -02	5.1003E -04
	90	2.1112E -08	2.8489E -02	1.9658E -05

(i) A direct numerical evaluation of the tangential derivatives of ϕ n in the integrands require, before all, the analytical computations of the Jacobian matrix of vector spherical harmonics and of the tangent and cotangent vectors (3.1)-(3.2) .

Acknowledgment Part of this work was carried out while the author was supported by DFG grant CRC 755 at the University of Goettingen. The author would thank Stéphanie Chaillat from the Laboratory POems, CNRS-INRIA Roquencourt-ENSTA ParisTech, for providing new numerical results [12] for comparison with our algorithm. The author is also grateful to the anonymous referees for their useful comments and suggestions.

Appendix

In this appendix, we give the analytical expressions of the solution to the time-harmonic scattering problems of a plane wave by an elastic sphere of radius R > 0 with a Dirichlet or a Neumann boundary condition.

Notations We denote by Y ℓ,j the spherical harmonics of order ℓ for j = -ℓ, . . . , ℓ, with ℓ, j ∈ N [30, Equation (2.4.78)]; by j ℓ the spherical Bessel function of order ℓ ∈ N and by h ℓ . We introduce the scalar function u

ℓ,j (κ s , x), k = 1, 3. These functions solve the Helmholtz equation for the wavenumber κ p and the Maxwell's equations for the wavenumber κ s [29, Theorems 9.14 and 9.16], respectively.

An incident plane wave defined by (5.1) admits the following series expansion (using [29, Equations (9.44), (9.45), (9.63)] and [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]Theorem 6.26])

with α

(1)

The vector spherical harmonics form a complete orthonormal system in the Hilbert space L 2 (S 2 ). Solving 3 × 3 linear systems associated to the boundary conditions (1.2) or (1.3) we obtain the series expansion of the scattered waves and then the one of the far-field pattern.

The Dirichlet boundary value problem For the rigid body problem we obtain, when the boundary of the elastic obstacle is a sphere of radius R > 0, the far-field is given by the following formula

where we have γ

and for ℓ ≥ 1 we have α

(1) ℓ

and for k = 1, 3

The Neumann boundary value problem For the cavity problem, when the boundary of the elastic obstacle is a sphere of radius R > 0, the far-field is given by the following formula

where we have γ

and for ℓ ≥ 1 we have

(1) ℓ and for k = 1, 3

and b