Ismaïl Badia
email: ismail.badia@doct.uliege.be

Boris Caudron
email: boris.caudron@thalesgroup.com

Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Christophe Geuzaine
email: cgeuzaine@uliege.be

A well-conditioned weak coupling of boundary element and high-order finite element methods for time-harmonic electromagnetic scattering by inhomogeneous objects

Keywords: Maxwell's equations, time-harmonic scattering, weak coupling, optimized domain decomposition method, high-order finite elements, boundary elements

The aim of this paper is to propose efficient weak coupling formulations between the boundary element method and the high-order finite element method for solving time-harmonic electromagnetic scattering problems. The approach is based on the use of a non-overlapping domain decomposition method involving optimal transmission operators. The associated transmission conditions are constructed through a localization process based on complex rational Padé approximants of the nonlocal Magnetic-to-Electric operators. Numerical results are presented to validate and analyze the new approach for both homogeneous and inhomogeneous scatterers.

Introduction

The time-harmonic Maxwell equations, which describe the propagation of electromagnetic waves, are essential for a wide range of technological applications as e.g. in antenna and radar design. Numerically solving three-dimensional time-harmonic electromagnetic scattering problems is known to be challenging since it leads to the discretization of a system of partial differential equations set in an unbounded domain. In addition, in the high frequency regime, i.e. when the wavelength λ is small compared to the characteristic size of the scatterer, very fine discretization grids are required to capture the highly oscillatory wave field solution, leading to a huge number of degrees of freedom and associated large size linear systems to solve.

Various approaches can be used to numerically solve time-harmonic electromagnetic scattering problems by a bounded dielectric and inhomogeneous object placed in a homogeneous domain which extends to infinity. Among the most widely used, we can mention the Finite Element Method (FEM) with an Absorbing/Artificial Boundary Condition (ABC) [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF] or a Perfectly Matched Layer (PML) [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Chew | A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF][START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF] to bound the computational domain. This method is well-suited for tackling complex geometrical configurations and heterogeneous media. However, using the FEM for wave problems is known to suffer from the so-called pollution effect [START_REF] Monk | A dispersion analysis of finite element methods for Maxwell's equations[END_REF] related to the accumulation of phase error over the mesh. From a practical point of view, this implies that the mesh size h must be adapted according to the wavenumber k = 2π/λ. A way to reduce the pollution error without refining the mesh is to increase the order of the shape functions involved in the FEM. In addition, solving the large size complex-valued and possibly indefinite linear system associated to the FEM approach, especially in the high-frequency regime, remains computationally challenging. To this end, specifically designed efficient and robust solvers are needed [START_REF] Alonso-Rodriguez | New nonoverlapping domain decomposition methods for the harmonic Maxwell system[END_REF][START_REF] Dolean | Optimized Schwarz methods for Maxwell's equations[END_REF][START_REF] Bouajaji | Optimized schwarz methods for the time-harmonic Maxwell equations with damping[END_REF][START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF][START_REF] Ernst | Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods[END_REF][START_REF] Peng | Non-conformal domain decomposition method with second-order transmission conditions for time-harmonic electromagnetics[END_REF][START_REF] Peng | One way domain decomposition method with second order transmission conditions for solving electromagnetic wave problems[END_REF][START_REF] Rawat | Nonoverlapping domain decomposition with second order transmission condition for the time-harmonic Maxwell's equations[END_REF][START_REF] Stupfel | High-order transmission conditions in a domain decomposition method for the time-harmonic Maxwell's equations in inhomogeneous media[END_REF]. To avoid the treatment of the exterior domain by the FEM, an alternative approach, called FEM-BEM coupling, consists in combining the FEM for the inhomogeneous part of the computational domain with a surface integral equation for the exterior domain, set on the boundary of the scatterer [START_REF] Nédélec | Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems[END_REF], and discretized by the Boundary Element Method (BEM). Although the discretization of an integral equation results in a complex-valued indefinite matrix which is fully populated and expensive to store and solve, advanced compression algorithms like the Fast Multipole Method (FMM) [START_REF] Liu | Fast Multipole Boundary Element Method: Theory and Applications in Engineering[END_REF] or Hierarchical Matrices (H-Matrices) [START_REF] Börm | Efficient Numerical Methods for Nonlocal Operators: H2-matrix Compression, Algorithms and Analysis[END_REF][START_REF] Hackbusch | Hierarchical Matrices: Algorithms and Analysis[END_REF] remedy these drawbacks. A FEM-BEM coupling therefore combines the advantages of both methods.

Standard FEM-BEM techniques consist in combining a surface integral equation for the exterior domain and a volume variational formulation for the interior domain within a single linear systeman approach often called strong FEM-BEM coupling (see e.g. [START_REF] De | Some formulations coupling finite element and integral equation methods for Helmholtz equation and electromagnetism[END_REF][START_REF] Hiptmair | Coupling of finite elements and boundary elements in electromagnetic scattering[END_REF][START_REF] Hiptmair | Stabilized FEM-BEM coupling for Maxwell transmission problems[END_REF][START_REF] Levillain | Couplage éléments finis-équations intégrales pour la résolution des équations de Maxwell en milieu hétérogène[END_REF][START_REF] Liu | A highly effective preconditioner for solving the finite element-boundary integral matrix equation of 3-D scattering[END_REF][START_REF] Sheng | On the formulation of hybrid finite element and boundary-integral methods for 3-D scattering[END_REF][START_REF] Yuan | Three-dimensional electromagnetic scattering from inhomogeneous objects by the hybrid moment and finite element method[END_REF]). Such direct formulations unfortunately have some major drawbacks that limit their application for solving highfrequency problems. First, strongly coupled methods yield very large size indefinite complex-valued linear systems defined by a matrix with has both sparse and dense parts. Such linear systems are generally difficult to solve and not adapted to compression algorithms for the nonlocal integral part. Moreover, it is not possible to combine two pre-existing efficient solvers, i.e. one FEM solver for the interior domain and one BEM solver for the exterior domain, to construct a global optimized solver for the coupled formulation. To remedy these drawbacks, alternative weak FEM-BEM couplings were recently introduced for the three-dimensional Helmholtz equation that arises in acoustics [START_REF] Bendali | A FETI-like domain decomposition method for coupling finite elements and boundary elements in large-size problems of acoustic scattering[END_REF][START_REF] Boubendir | Coupling of a non-overlapping domain decomposition method for a nodal finite element method with a boundary element method[END_REF][START_REF] Caudron | Optimized weak coupling of boundary element and finite element methods for acoustic scattering[END_REF] and for the time-harmonic Maxwell's equations [START_REF] Stupfel | High-order transmission conditions in a domain decomposition method for the time-harmonic Maxwell's equations in inhomogeneous media[END_REF]. Non-conformal weak FEM-BEM couplings have also been proposed [START_REF] Peng | Non-conformal domain decomposition method with second-order transmission conditions for time-harmonic electromagnetics[END_REF][START_REF] Vouvakis | A domain decomposition approach for nonconformal couplings between finite and boundary elements for unbounded electromagnetic problems in R 3[END_REF] (e.g. independently generated meshes for different parts of the problem under analysis) but may introduce additional approximation errors. The weak FEM-BEM approach can be interpreted as an optimized Schwarz domain decomposition method (DDM) iterating between the bounded scatterer and the exterior complementary domain until a global convergence criterion is reached. One of the strengths of the weak FEM-BEM approach is that using two pre-existing optimized solvers, one for the BEM and one for the FEM, is direct. Therefore, in the context of an industrial collaboration, using two independently developed solvers could be a considerable advantage. Indeed, for each interior/exterior electromagnetic subproblem, different partners (e.g. aircraft/antenna designers) can exchange minimal information at the transmitting boundary between the two subdomains to solve the global electromagnetic problem without sharing their simulation codes and the details of the models.

To design an efficient weak FEM-BEM coupling, it is crucial to consider a fast converging Schwarz DDM algorithm. For wave propagation problems, and most particularly in the highfrequency regime, the key computational ingredient in the implementation of an efficient and robust converging DDM is to design well-suited transmission conditions between the interfaces of the subdomains. In electromagnetism, the optimal transmission condition is defined through the Magnetic-to-Electric (MtE) map that links the magnetic and the electric surface currents at the interface of the subdomains [START_REF] Caudron | Couplages FEM-BEM faibles et optimisés pour des problémes de diffraction harmoniques en acoustique et en électromagnétisme[END_REF]. From a computational point of view, the numerical evaluation of this operator is nevertheless costly since the MtE operator is a nonlocal pseudodifferential operator. This is the reason why Schwarz DDM with optimized local transmitting boundary conditions were proposed in the literature [START_REF] Alonso-Rodriguez | New nonoverlapping domain decomposition methods for the harmonic Maxwell system[END_REF][START_REF] Després | A domain decomposition method for the harmonic Maxwell equations[END_REF][START_REF] Dolean | Optimized Schwarz methods for Maxwell's equations[END_REF][START_REF] Bouajaji | Optimized schwarz methods for the time-harmonic Maxwell equations with damping[END_REF][START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF][START_REF] Rawat | Nonoverlapping domain decomposition with second order transmission condition for the time-harmonic Maxwell's equations[END_REF]. In the context of the weak FEM-BEM coupling, it was shown in [START_REF] Caudron | Couplages FEM-BEM faibles et optimisés pour des problémes de diffraction harmoniques en acoustique et en électromagnétisme[END_REF][START_REF] Caudron | Optimized weak coupling of boundary element and finite element methods for acoustic scattering[END_REF] that a fast convergence can be expected when using the Padé expansion of a square-root operator which approximates the nonlocal Dirichlet-to-Neumann (DtN) operator, where the DtN map plays the role of the MtE operator but for acoustics. The aim of the paper is to present the extension and some improvements of the weak FEM-BEM coupling method to the electromagnetic scattering problem by an inhomogeneous scatterer by using high-order finite element methods. This new weak FEM-BEM coupling formulation is shown to exhibit a convergence rate that is only slightly dependent on the wavenumber k, the mesh refinement h and the contrast between the two subdomains. Note that for homogeneous obstacles one can use a weak BEM-BEM coupling formulation.

The article is organized as follows. Section 2 describes the electromagnetic transmissionscattering problem. In Section 3, we present the weak coupling approach. Section 4 is devoted to the construction of approximate nonlocal MtE surface operators. Section 5 provides the formulations for solving the weak coupling and associated subproblems. Section 6 details the finite element implementation of the weak coupling. Numerical results for three-dimensional electromagnetic scattering problems by homogeneous and inhomogeneous obstacles are presented in Section 7. We finally conclude in Section 8.

The electromagnetic transmission-scattering problem

Let Ω -be a bounded scatterer in R 3 with smooth closed boundary Γ := ∂Ω -. The exterior unbounded domain of propagation is defined by Ω + := R 3 \ Ω -. We consider an incident electromagnetic plane wave (E inc ;H inc ) propagating in the unbounded domain Ω + and satisfying the exterior time-harmonic Maxwell's equations

curl E inc -ιk + Z + H inc = 0 curl H inc + ιk + Z -1 + E inc = 0 in Ω + .
The object generates an outgoing scattered wave (E sc ; H sc ). Let us define (E + ; H + ) = (E inc ; H inc )+ (E sc ; H sc ) as the total electromagnetic field in the exterior domain Ω + . Then, our goal is to solve the following boundary-value problem: find the wave field E ∈ H(curl , Ω -) ∪ H loc (curl , Ω + ) 1 solution of the following electromagnetic transmission-scattering problem

curl 1 k -Z - curl E -k -Z -1 -E = 0, in Ω -, (2.1a
)

curl curl E -k 2 + E = 0, in Ω + , (2.1b) 1 ιk + curl E sc × x x -E sc = O r→+∞ (r -2 ), (2.1c)
γ - t E = γ + t E , on Γ, (2.1d) γ - t 1 k -Z - curl E = 1 k + Z + γ + t (curl E ) , on Γ. ( 2.1e)
In the above equations,

k ± := ω √ ε ± µ ± and Z ± := µ ±
ε ± denote respectively the wavenumbers and the impedances associated with Ω ± , where ω is the angular frequency, and ε ± and µ ± stand for the permittivities and permeabilities in Ω ± . In practice, the exterior domain is made of air so that Z + corresponds to the free space impedance Z 0 . To include the case of heterogeneous dielectric materials, k -and Z -can be space-dependent functions of x ∈ R 3 . We introduce a contrast parameter δ ∈ R such as, in case of homogeneous scatterers and µ + = µ -, k -= δk + and Z -= Z + δ . The curl operator of a vector field u ∈ R 3 is denoted by curl u. The unit imaginary number is ι = √ -1. The notation u × v designates the cross product and u • v the inner product between two complex valued vector fields u and v in C 3 , where v is the complex conjugate of v. The associated norm is x := √ x • x. Equation (2.1c) is the Silver-Müller radiation condition at infinity with x = rσ ∈ R 3 , r := x the radial coordinate and σ the unit radial vector. Finally, the transmission conditions (2.1d) and (2.1e) respectively express the tangential continuity of the total electric and magnetic fields through Γ. These transmission conditions involve the tangential trace operators γ ± t defined by

γ ± t : H(curl , Ω -) ∪ H loc (curl , Ω + ) → H -1 2 t (div Γ , Γ) u → u ±|Γ × n ,
with n the outward-pointing unit normal vector of Ω -on Γ and H

-1 2 t (div Γ , Γ) = {v ∈ H -1 2 t (Γ) : div Γ u ∈ H -1 2 (Γ)}
, where the surface divergence operator div Γ [START_REF] Nédélec | Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems[END_REF] has to be taken in the distributional sense. The space H -1 2 t (Γ) designates the dual of the Hilbert space H

1 2 t (Γ) = {u ∈ H 1 2 (Γ) : 1 H(curl , Ω-) = {v ∈ L 2 (Ω-) : curl v ∈ L 2 (Ω-)} ; H loc (curl , Ω+) = {v ∈ H(curl , B \ Ω-) :
for every open ball B containing Ω-in its interior}. u • n = 0}. Provided that k -and Z -are piecewise smooth and bounded away from zero everywhere in Ω -, the transmission-scattering problem (2.1) has a unique solution [START_REF] Hiptmair | Coupling of finite elements and boundary elements in electromagnetic scattering[END_REF]. An illustration of the problem is given in Fig. 1.

The weak coupling approach

A classical way to solve (2.1) is to couple a surface integral formulation for the exterior problem and a variational volume formulation for the interior domain within a single equation, leading to a strong FEM-BEM coupling. As noted in the introduction, these methods have however several drawbacks. We therefore rather consider a weak coupling approach. This method is essentially a non-overlapping Schwarz DDM based on the following decomposition of R 3 : R 3 = Ω -∪ Ω + ∪ Γ. The Schwarz DDM is solved in practice via an iterative method such as GMRES [START_REF] Saad | Iterative Methods for Sparse Linear systems[END_REF][START_REF] Saad | GMRES: A Generalized Minimal Residual algorithm for solving nonsymmetric linear systems[END_REF]. The key point to get a fast converging DDM is to design optimized transmission conditions at the interface Γ. For the time-harmonic Maxwell equations, Després first proposed in [START_REF] Després | A domain decomposition method for the harmonic Maxwell equations[END_REF] to use a simple impedance boundary operator to obtain a converging algorithm. Several other local transmission conditions were later derived (see e.g. [START_REF] Alonso-Rodriguez | New nonoverlapping domain decomposition methods for the harmonic Maxwell system[END_REF][START_REF] Després | A domain decomposition method for the harmonic Maxwell equations[END_REF][START_REF] Dolean | Optimized Schwarz methods for Maxwell's equations[END_REF][START_REF] Bouajaji | Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell's equations[END_REF][START_REF] Bouajaji | Optimized schwarz methods for the time-harmonic Maxwell equations with damping[END_REF][START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF][START_REF] Peng | One way domain decomposition method with second order transmission conditions for solving electromagnetic wave problems[END_REF][START_REF] Rawat | Nonoverlapping domain decomposition with second order transmission condition for the time-harmonic Maxwell's equations[END_REF]). Similarly, the well-conditioned weak coupling algorithm requires the use of adapted interface conditions. In this work, we will use nonlocal MtE maps and their localization through Padé approximants [START_REF] Caudron | Couplages FEM-BEM faibles et optimisés pour des problémes de diffraction harmoniques en acoustique et en électromagnétisme[END_REF][START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] to achieve fast convergence.

Reformulating the problem

The first step of the method consists in reformulating the transmission conditions (2.1d)-(2.1e). Let us introduce the following two transmission operators

T -: H -1 2 t (div Γ , Γ) → H -1 2 t (div Γ , Γ), T + : H -1 2 t (div Γ , Γ) → H -1 2 t (div Γ , Γ). (3.1)
Since we work with total fields, T ± are affine operators. We assume that T -+T + is injective so that the transmission conditions (2.1d)-(2.1e) are equivalent to the following transmission conditions

γ - t H + T -γ - t E = γ + t H + T -γ + t E γ + t H -T + γ + t E = γ - t H -T + γ - t E . (3.2)
DDMs are usually constructed by assuming that T + = T -(see e.g. [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF]). Here, however, we consider that they can be different for more generality and flexibility. Then, we introduce the following trace operators

∀(E ; H) ∈ [H(curl , Ω -) ∪ H loc (curl , Ω + )] 2 , B ± (E ; H) = γ ± t H ∓ T ± γ ± t E,
from which we define the unknowns of the weak coupling

g -= B -(E; H), g + = B + (E; H). (3.3)
Note that by definition of the unknowns of the weak coupling (3.3), it follows that g ± = 0. Finally, we still need to introduce the affine resolution operators R -and R + as: ∀g ∈ H

-1 2 t (div Γ , Γ), R ± g = γ ± t E ± ,
   curl E --ιk -Z -H - = 0 in Ω - curl H -+ ιk -Z -1 -E - = 0 in Ω - B -(E -; H -) = g on Γ , (3.4)
and

         curl E + -ιk + Z + H + = 0 in Ω + curl H + + ιk + Z -1 + E + = 0 in Ω + Z + H sc × x ||x|| -E sc = O r→+∞ (r -2 ) B + (E + ; H + ) = g on Γ . ( 3.5)
The subproblems (3.4) and (3.5) are associated with the DDM for the weak coupling. Their wellposedness depends on the choice of the transmission operators. We assume here that both problems (3.4) and (3.5) are well-posed so that the resolution operators are correctly defined. We can now introduce the weak coupling formulation of the transmission-scattering problem (2.1).

Proposition 3.1. The electromagnetic field (E; H) is the solution of the transmission-scattering problem (2.1) if and only if g -and g + are solutions to

(Id -S π ) g - g + = 0 0 , (3.6)
the operator Id -S π being defined by

Id -S π = Id -S + -S -Id , with S ± = Id ± (T -+ T + )R ± .
Proof. We start by rewriting the transmission conditions (3.2) as

B -(E, H) = B + (E, H) + (T -+ T + )γ + t E, B + (E, H) = B -(E, H) -(T -+ T + )γ - t E.
Since the field (E; H) is solution of (3.4) and (3.5), we deduce that

R ± g = γ ± t E,
and obtain (3.6). Finally to prove the equivalence, it is sufficient to check that the transmission conditions (3.2) are satisfied, which can be easily proved.

In a domain decomposition setting, let us remark that an iterative Jacobi algorithm [START_REF] Saad | Iterative Methods for Sparse Linear systems[END_REF] for solving (3.6) involves at iteration n + 1 the resolution of the following subproblems

   curl E n+1 - -ιk -Z -H n+1 - = 0 in Ω - curl H n+1 - + ιk -Z -1 -E n+1 - = 0 in Ω - B -(E n+1 -; H n+1 -) = g n - , (3.7)
and

         curl E n+1 + -ιk + Z + H n+1 + = 0 in Ω + curl H n+1 + + ιk + Z -1 + E n+1 + = 0 in Ω + Z + H n+1 sc × x ||x|| -E n+1 sc = O r→+∞ (r -2 ) B + (E n+1 + ; H n+1 + ) = g n + , (3.8)
and finally the updating of the interface unknowns according to

g n+1 - g n+1 + = S π g n - g n + . (3.9)
This is equivalent to the following classical Schwarz algorithm

             curl E n+1 -ιk ± Z ± H n+1 = 0 in Ω ± curl H n+1 + ιk ± Z -1 ± E n+1 = 0 in Ω ± Z + H n+1 sc × x ||x|| -E n+1 sc = O r→+∞ (r -2 ) γ - t H n+1 + T -γ - t E n+1 = γ + t H n + T -γ + t E n γ + t H n+1 -T + γ + t E n+1 = γ - t H n -T + γ - t E n .
In (3.7) and (3.8), we can interpret g n -and g n + as some exchanged information from Ω + to Ω -and from Ω -to Ω + .

Optimal transmission operators

The weak coupling (3.6) is not suitable for a direct resolution. Indeed, the operator Id -S π is defined implicitly from the resolution operators R + and R -, which correspond, at the discrete level, to solving (3.4) and (3.5), respectively. In this work, the final weak coupling linear system is solved via an iterative Krylov subspace method (GMRES) [START_REF] Saad | Iterative Methods for Sparse Linear systems[END_REF][START_REF] Saad | GMRES: A Generalized Minimal Residual algorithm for solving nonsymmetric linear systems[END_REF].

The convergence of the weak coupling is fundamentally related to the choice of the transmission operators T ± . To this end, let us introduce the affine nonlocal MtE operators Λ ±,k ± ,Z ± defined by

Λ ±,k ± ,Z ± : H -1 2 t (div Γ , Γ) → H -1 2 t (div Γ , Γ) γ t ± E → γ t ± H .
By rewriting (3.3), the MtE, transmission and resolution operators are linked by

(Λ -,k -,Z -+ T -)R -= Id, (Λ +,k + ,Z + -T + )R + = Id.
Let us now consider the following important proposition.

Proposition 3.2. The optimal choice of the transmission operators for the weak coupling (3.6) is given as

T -= -Λ +,k + ,Z + ⇒ S + = 0 , T + = Λ -,k -,Z -⇒ S -= 0. (3.10)
A direct consequence of this result is that S π is null when using the MtE operators as transmission operators, leading to a trivial resolution of the weak coupling (3.6). In addition, if only one of the transmission operators T ± is optimal, i.e. one of the two conditions (3.10) holds, then the linear part of S π is nilpotent and its spectrum is therefore reduced to {0}. Therefore, the linear part of Id -S π is characterized by an eigenvalue clustering on (1, 0) in the complex plane, leading to a fast convergence of the GMRES used to solve the weak coupling. The MtE operators are a priori nonlocal operators that lead to high computational costs. To get a reduced evaluation cost for the weak coupling procedure, it is highly desirable to rather consider localized approximations of the MtE operators. Such accurate representations are available in the literature for the exterior MtE operator Λ +,k + ,Z + (see Section 4). This can be understood by the fact that the solution to the exterior problem is spatially localized at least for convex domains, in particular in the high-frequency regime, where constructive methods are then possible. Therefore, it is reasonable to get access to an accurate computation of T -= -Λ +,k + ,Z + , which then leads to an almost nilpotent linear part of S π . Building an accurate approximation of T + = Λ -,k -,Z -is more difficult, and can even be impossible by localized operators due to the presence of multiply bouncing internal waves in Ω - corresponding to complex nonlocal interactions. Since we can get an accurate representation of T - and almost a nilpotent linear part of S π , a rough approximation of Λ -,k -,Z -to get T + , e.g. by a constant, may seem to be sufficient. However, this can sometimes not be the case as it was observed in [START_REF] Caudron | Optimized weak coupling of boundary element and finite element methods for acoustic scattering[END_REF] for the acoustic case and here for the electromagnetic case (see Section 7) where an increase of GMRES iterations can be observed at some specific frequencies, while leading nevertheless to a converging weak coupling procedure with a correct solution. This drawback is accentuated when the contrast between the exterior and interior domains is high. An alternative choice that appears first as less natural for a two domain decomposition is to rather extend in a fictive way the exterior problem into Ω -by setting T + = -Λ +,k + ,Z + = T -. This leads to improved convergence at these specific frequencies, even with a rough approximation of T + by a constant (and then T + = T - when localized). We will see in Section 7 that this improvement is important most particularly for the high contrast case.

Finally, let us remark that an alternative choice of the operators T ± could have been made to optimize the resolution of the subproblems (3.4) and (3.5), e.g. to get well-conditioned integral equation formulations [START_REF] Caudron | Couplages FEM-BEM faibles et optimisés pour des problémes de diffraction harmoniques en acoustique et en électromagnétisme[END_REF]. These optimal operators are different from those of the weak coupling. Since our main objective is to precondition the weak coupling (3.6), we only consider here the adapted transmission operators (3.10).

Approximations of the MtE maps

The MtE operators Λ ±,k,Z are usually not available for a general shape Γ. Let us note that for a magnetic current m, we can write m = m + m inc where m inc is the contribution of the incident field (m inc = 0 if there is none) and m = mm inc . In the following, we propose to decompose the affine operators

Λ ±,k,Z such that ∀m, Λ ±,k,Z (m) = Λ ±,k,Z ( m) + Λ inc ±,k,Z
where Λ ±,k,Z are linear operators and Λ inc ±,k,Z = Λ ±,k,Z (m inc ) (for the contribution of the incident field).

Nonlocal approximations of the linear part of the MtE operators

Let us assume that Γ is a general smooth convex surface. Based on microlocal analysis [START_REF] Bouajaji | Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell's equations[END_REF], we can build the following nonlocal surface approximations

Λ sq ±,k,Z of Λ ±,k,Z Λ sq ±,k,Z = ∓ 1 Z (Id + ∇ Γ ( 1 k 2 div Γ ) -curl Γ ( 1 k 2 curl Γ )) -1 2 (Id -curl Γ ( 1 k 2 curl Γ ))(Id × n), (4.1)
with k = k + ι . The accuracy of (4.1) used as an OSRC [START_REF] Bouajaji | Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell's equations[END_REF] depends crucially on the choice of the regularization parameter . A well-suited value is = 0.39k

1 3 H 2 3
, with H the local mean curvature. Since the operators Λ sq ±,k ± ,Z ± are used as preconditioners for the weak coupling procedure, we can simplify this expression by replacing H by the inverse of the radius of the circumscribed sphere to Ω -. The operators ∇ Γ , curl Γ and curl Γ designate the surface gradient, the vector and scalar surface curl operators [START_REF] Nédélec | Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems[END_REF], respectively. If k and Z are space-dependent functions, we will formally extend (4.1) which can be justified thanks to the theory of pseudodifferential operators and the associated microlocal analysis techniques (see e.g. some examples in [START_REF] Ph | Construction and numerical assessment of local absorbing boundary conditions for heterogeneous time-harmonic acoustic problems[END_REF] for scalar wave equation). It is worth noting that the theory of pseudodifferential operators is developed by assuming that the surface is smooth. By weakening the regularity, for example near a corner, the approximations (4.1) may lose some accuracy since they are not rigorously valid. However, again because they numerically act as preconditioners for the weak coupling algorithms, this has a limited impact in the global convergence of the iterative scheme (e.g. GMRES). It would however be very interesting to improve the convergence of the current method by carefully including the corner effects [START_REF] Modave | Corner treatments for high-order local absorbing boundary conditions in high-frequency acoustic scattering[END_REF].

In practice, Λ sq +,k + ,Z + provides a suitable approximation of Λ +,k + ,Z + , and then the linear part of T -for the weak coupling, in particular for large wavenumbers k + . However, since Λ -,k -,Z -is usually a strongly nonlocal operator, its approximation by Λ sq -,k -,Zis relatively limited. For k - and Z -constant, a possible alternative would be to use an integral approximation (see e.g. Section 5.3.3, pages 133-134 in [START_REF] Caudron | Couplages FEM-BEM faibles et optimisés pour des problémes de diffraction harmoniques en acoustique et en électromagnétisme[END_REF]). However, this choice may appear as costly in practice. Since we are already considering a suitable approximation of T -, as mentioned before, a rough approximation of T + is however sufficient because the linear part of S π is almost nilpotent. We will see in Section 7 that an approximation of the linear part of T + by a brute force approximation of -Λ sq +,k + ,Z + based on a constant leads to very satisfactory convergence results.

Localization of the linear part of the MtE operators

The approximations (4.1) involve nonlocal pseudodifferential operators that are computationally expensive to evaluate directly and lead to full discrete matrix representations (when they can be computed). A practical approach is to further use localized approximations of the inverse squareroot symbol. The simplest approximation of Λ sq +,k,Z and Λ sq -,k,Z is to use a Taylor expansion in 1/k truncated at the zeroth-order term following

Λ 0 ±,k,Z = ∓ 1 Z (Id × n). (4.2)
More accurate local approximations of Λ sq ±,k,Z , leading to Padé-localized square-root approximations that we denote by Λ sq,Np,θp ±,k,Z , were proposed in [START_REF] Bouajaji | Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell's equations[END_REF]. To this end, let us introduce the rational Padé approximation of order N p of the square-root function [START_REF] Milinazzo | Rational square-root approximations for parabolic equation algorithms[END_REF], with θ p -rotation of the branch-cut:

∀z ∈ C, R(z) > -1, ∀j ∈ {1; ...; N p }, z = -B -1 j , (1 + z) 1 2 ≈ e ι θp 2 R Np ((1 + z)e -ιθp -1) = C 0 + Np =1 A z 1 + B z , (4.3)
where R Np is the usual real-valued Padé approximation of order N p

R Np = 1 + Np =1 a z 1 + b z with a = 2 2N p + 1 sin 2 π 2N p + 1 and b = cos 2 π 2N p + 1 .
The angle of rotation θ p is a free parameter that is fixed for the numerical simulations and

C 0 = e ι θp 2 R Np (e -ιθp -1), A = e -ι θp 2 a (1 + b (e -ιθ -1)) 2 , B = e -ιθp b 1 + b (e -ιθ -1)
.

Formally, considering that

J = ∇ Γ ( 1 k 2 div Γ ) -curl Γ ( 1 k 2 curl Γ ), the Padé approximation (4.
3) leads to the following approximate linear operators

Λ sq,Np,θp ±,k,Z of Λ sq ±,k,Z Λ sq,Np,θp ±,k,Z = ∓ 1 Z   C 0 Id + Np =1 A J (Id + B J ) -1   -1 Id -curl Γ ( 1 k 2 curl Γ ) (Id × n). (4.4)
Since the operators (4.4) are expressed thanks to local and inverse of local surface operators, the implementation of Λ sq,Np,θp +,k + ,Z + can be realized locally by introducing some auxiliary surface fields (see Section 5.2). This is an important point since only complex-valued sparse matrices appear at the discrete level in the weak coupling procedure when using -Λ sq,Np,θp

+,k + ,Z + to approximate -Λ +,k + ,Z + .
This results in an additional low memory storage as well as an efficient implementation when applying this approximate linear operator.

Since T -is accurately evaluated, we can consider a simple approximation for T + . In particular, it was shown in [START_REF] Caudron | Couplages FEM-BEM faibles et optimisés pour des problémes de diffraction harmoniques en acoustique et en électromagnétisme[END_REF] that, in the high-frequency regime, the convergence properties are very similar when using Padé-type or low-order approximations of the linear part of T + . Therefore, we restrict our study to the zeroth-order approximations (4.2) for the linear part of T + . In addition, choosing Λ 0 -,k -,Zappears as natural. However, it was shown in the acoustic case [START_REF] Caudron | Couplages FEM-BEM faibles et optimisés pour des problémes de diffraction harmoniques en acoustique et en électromagnétisme[END_REF] that the number of GMRES iterations can still be high at some frequencies, which is amplified for high contrast situations. For this reason, we also propose to consider the alternative choice Λ 0 -,k + ,Z + which provides an improved convergence rate as validated in Section 7.

Formulations for solving the weak coupling and associated subproblems

We now present the weak coupling formulations. For the interior problem (3.4), we use a variational volume formulation which is more adapted to simulate non homogeneous materials when discretized by high-order finite elements. A boundary element method is considered for the exterior problem (3.5). To this end, we first recall some basic results about the integral representation for Maxwell's equations.

Basics on integral representation

Let us consider a wavenumber k ∈ R * + and an impedance Z ∈ R * + . We denote by G k the three-dimensional Green's function defined by

∀x ∈ R 3 \ {0}, G k (x) = e ιk x 4π
x .

The electric and magnetic field potential operators are respectively given by ∀j ∈ H -1

and

K k j(x) = -curl Γ G k (x -y)j(y) dy.
Taking the traces of T k j and K k j on Γ and introducing the electric boundary operator: ∀j ∈ H -1 2 t (div Γ , Γ), and ∀x ∈ Γ,

T k j(x) = Γ G k (x -y)j(y) dy × n(x) + 1 k 2 ∇ Γ G k (x -y)div Γ j(y) dy × n(x)
and the magnetic boundary operator: ∀j ∈ H -1 2 t (div Γ , Γ) and ∀x ∈ Γ,

K k j(x) = Γ ∇ y G k (x -y) × j(y) dy × n(x),
we obtain the following results.

Proposition 5.1. We have the trace relations:

γ - t T k = T k , γ + t T k = T k , γ - t K k = - 1 2 Id + K k , γ + t K k = 1 2 Id + K k .
To derive well-posed integral equation formulations, we need to introduce the integral representation for Maxwell's equations [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF][START_REF] Nédélec | Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems[END_REF].

Proposition 5.2. Introducing the so-called jump relations across Γ

[γ t E] := γ - t E -γ + t E sc and [γ t H] := γ - t H -γ + t H sc
, and assuming that Ω -is a homogeneous domain, we then have the following integral representation

E = E inc + ιkZT k [γ t H] -K k [γ t E] H = H inc -K k [γ t H] -ιkZ -1 T k [γ t E] in Ω -∪ Ω + , (5.1
)

with E inc = 0 H inc = 0 in Ω -.
Once the jump relations are known or estimated on Γ, the representation formulae (5.1) can be used to find the solutions in Ω + and in a homogeneous interior domain Ω -.

Formulations for solving the weak coupling and the exterior/interior subproblems

Several possible combinations of functional spaces can be considered for g ± :

(g -, g + ) ∈ H -1 2 t (div Γ , Γ) × H -1 2 t (rot Γ , Γ), (g -, g + ) ∈ H -1 2 t (rot Γ , Γ) × H -1 2 t (rot Γ , Γ) and (g -, g + ) ∈ H -1 2 t (rot Γ , Γ) × H -1 2 t (div Γ , Γ).
From the functional analysis point of view, the best choice is not obvious. During the numerical simulations in Section 7, it however appears that choosing g ± ∈ H -1 2 t (div Γ , Γ) leads to the best convergence properties. Let us mention that alternative approaches such as the Helmholtz decomposition [START_REF] Caudron | Couplages FEM-BEM faibles et optimisés pour des problémes de diffraction harmoniques en acoustique et en électromagnétisme[END_REF] could also be used. The spaces H -1 2 t (rot Γ , Γ) and H -1 2 t (div Γ , Γ) are mutually adjoint with respect to the scalar product L 2 t (Γ) [START_REF] Nédélec | Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems[END_REF]. In preparation for the discrete formulation in Section 6, we will however formally identify, in the following proposition, the H

-1 2 t (div Γ , Γ) space with the discrete H -1 2 t (div Γ , Γ) space which is L 2 t (Γ) with itself. Proposition 5.3. Let us consider some test-functions g -∈ H -1 2 t (div Γ , Γ) and g + ∈ H -1 2
t (div Γ , Γ). Then, the variational formulation of the weak coupling (3.6) is given by

     Γ g -• g -dΓ - Γ g + • g -dΓ - Γ (T -+ T + )R + g + • g -dΓ = 0 Γ g + • g + dΓ - Γ g -• g + dΓ + Γ (T -+ T + )R -g -• g + dΓ = 0 .
(5.2)

The above weak formulation involves the resolution operators R ± associated with the subproblems (3.4) and (3.5). Since the physical parameters of the exterior problem (3.5) are supposed to be constant, the use of surface integral equations is particularly well-adapted. Because the exterior subproblem involves an impedance boundary condition, we consider some specific integral equation formulations. The chosen formulation corresponds to a generalization of the one proposed in [START_REF] Bachelot | Time dependent integral method for Maxwell's system with impedance boundary condition[END_REF].

Proposition 5.4. Let g ∈ H -1 2 t (div Γ , Γ) and (E; H) ∈ H loc (curl , Ω + ) 2 such that          curl E -ιk + Z + H = 0 in Ω + curl H + ιk + Z -1 + E = 0 in Ω + Z + H sc × x ||x|| -E sc = O r→+∞ (r -2 ) γ + t H -T + γ + t E = g .
We then have:

1 2 T -1 + γ + t H -K k + γ + t E + ιk + Z + T k + γ + t H = 1 2 T -1 + g + γ + t E inc 1 2 T + γ + t E -K k + γ + t H -ιk + Z -1 + T k + γ + t E = -1 2 g + γ + t H inc . (5.3)
Conversely, if g ∈ H -1 2 t (div Γ , Γ), m ∈ H -1 2 t (div Γ , Γ), j ∈ H -1 2 t (div Γ , Γ) satisfy 1 2 T -1 + j -K k + m + ιk + Z + T k + j = 1 2 T -1 + g + γ + t E inc 1 2 T + m -K k + j -ιk + Z -1 + T k + m = -1 2 g + γ + t H inc
, then the field (E; H) ∈ H loc (curl , Ω + ) 2 defined by

E = E inc -ιk + Z + T k + j + K k + m H = H inc + K k + j + ιk + Z + -1 T k + m in Ω +
is solution of the exterior problem (3.5).

Proof. Let (E; H) ∈ H loc (curl , Ω + ) 2 be the solution to the exterior problem (3.5). Extending (E; H) by (-E inc ;-H inc ) in Ω -, one gets the following integral representation formulae

E = E inc -ιk + Z + T k + γ + t H + K k + γ + t E H = H inc + K k + γ + t H + ιk + Z + -1 T k + γ + t E in Ω + .
Taking the tangential traces (5.1) and using the boundary condition, one then obtains (5.3). To prove the equivalence, it is sufficient to remark that the tangential traces of (E; H) are given by

γ + t E = γ + t E inc -ιk + Z + T k + j + K k + m + 1 2 m γ + t H = γ + t H inc + 1 2 j + K k + j + ιk + Z -1 + T k + m .
A calculation then shows that the boundary condition is fulfilled.

The interior subproblem is solved using a volume variational formulation. t (div Γ , Γ) and itself, the weak volume variational formulation for E associated with the interior problem (3.4) is ∀v ∈ H(curl, Ω -),

Ω - curl E • curl v dΩ -- Ω - k 2 -E • v dΩ - -ιk -Z -γ - t v; T -γ - t E = -ιk -Z -γ - t v; g .
(5.4)

The variational formulations of the linear part of the transmission operators based on (4.2) and (4.4) are presented in the following proposition.

Proposition 5.6. For all test functions v ∈ H -1 2 t (div Γ , Γ), the weak formulations associated with the linear part of the zeroth-order Taylor transmission conditions and the Padé-localized square-root transmission conditions are respectively defined for j ∈ H

-1 2 t (div Γ , Γ) by Γ Λ 0 ±,k,Z (j) • v dΓ = Γ ∓ 1 Z j × n • v dΓ (5.5) and Γ Λ sq,Np,θp ±,k,Z (j) • v dΓ = Γ ∓ 1 Z r k • v dΓ. (5.6)
The function

r k ∈ H -1 2 t (curl Γ , Γ
) is obtained by solving the following system with auxiliary variables: find (r k , (Ψ ) ∈{1;...;Np} , (φ ) ∈{1;...;Np} ) ∈ W such that, ∀(r , (Ψ ) ∈{1;...;Np} , (φ ) ∈{1;...;Np} ) ∈ W ,

                                 Γ C 0 r k • r dΓ - Γ j × n • r dΓ + Γ 1 k 2 curl Γ (j × n) • curl Γ (r ) dΓ + Np =1 A Γ ∇ Γ φ • r dΓ - Γ 1 k 2 curl Γ Ψ • curl Γ (r ) dΓ = 0 Γ Ψ • Ψ dΓ + B Γ ∇ Γ φ • Ψ dΓ - Γ 1 k 2 curl Γ Ψ • curl Γ (Ψ ) dΓ - Γ r k • Ψ dΓ = 0 ∀ ∈ {1; ...; N p } Γ φ • φ dΓ + Γ 1 k 2 Ψ • ∇ Γ φ dΓ = 0 ∀ ∈ {1; ...; N p } ,
where we set

W := H -1 2 t (curl Γ , Γ) × [H -1 2 t (curl Γ , Γ)] Np × [H 1 (Γ)] Np .

Finite element implementation of the weak coupling

Let us now explain the implementation of the electromagnetic weak FEM-BEM coupling. Its discretization is significantly more complex than the acoustic weak FEM-BEM coupling [START_REF] Caudron | Optimized weak coupling of boundary element and finite element methods for acoustic scattering[END_REF] due to the fact that the unknowns g ± are vectorial. More precisely, the main difficulty when discretizing the electromagnetic weak FEM-BEM coupling consists in choosing the representation of g ± at the discrete level.

Boundary/Finite element discretization

In this paper, the discretization of the surface integral equations is based on a Galerkin boundary element method (BEM). Other approaches can be used as e.g. the collocation and Nyström techniques [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF]. Unlike FEM, BEM enjoys a low numerical dispersion because the fundamental solution of the Maxwell's equations is explicitly taken into account in the integral formulations. Therefore, low-order approximation polynomials (typically linear elements) can be used. Conversely, it is well-known that dispersion/pollution errors appear in the numerical solution of FEM when increasing the frequency [START_REF] Monk | A dispersion analysis of finite element methods for Maxwell's equations[END_REF]. As a consequence, the mesh size h strongly depends on the wavenumber k. Using this strategy leads to considering very fine meshes then limiting the use of the standard low-order polynomial basis functions. An alternative way to reduce the error without refining the mesh is to increase the order of the shape functions, where the number of unknowns can still be controlled while keeping an accurate solution. Here, we use such a strategy based on high-order hierarchical finite elements [START_REF] Solín | Higher-order Finite Element Methods[END_REF]. Among hierarchic shape functions, one of the most popular choices is Lobatto shape functions. They are split into four different types of shape functions: vertex, edge, face and bubble functions. There are several benefits in using the hierarchical shape functions. Indeed, it was shown [START_REF] Solín | Higher-order Finite Element Methods[END_REF] that they provide good performances while being wellconditioned. Furthermore, their hierarchical nature implies that the shape functions at order p + 1 are obtained from the shape functions at order p by only adding new shape functions. Thus, one can very simply couple high-and low-order finite elements as well as apply static condensation for bubble shape functions, which vanish on the boundary of each element.

Let us assume that Ω -is approximated by a computational domain Ω h using n t tetrahedral finite elements with n f faces, n e edges and n v vertices. Then, Γ h , which is the boundary of Ω h , is composed of n ∂ t triangles, n ∂ e edges and n ∂ v vertices. In order to discretize the weak-coupling, we introduce some finite dimensional approximation subspaces such that Z h,p ⊂ H

-1 2 t (curl Γ , Γ), V h,p ⊂ H -1 2 t (div Γ , Γ), W h,p ⊂ H(curl, Ω -) and Q h ⊂ H 1 (Γ) with Z h,p := span{(z p j ) 1≤j≤N p z }, where N p z =    (p + 1)n ∂ e p < 2 (p + 1)n ∂ e + 3(p -1)n ∂ t 2 ≤ p < 3 (p + 1)n ∂ e + (p -1)(p + 1)n ∂ t 3 ≤ p , V h,p := span{(e p j ) 1≤j≤N p z } = span{(z p j × n) 1≤j≤N p z }, W h,p := span{(w j ) 1≤j≤N p w },
where

N p w =        (p + 1)n e p < 2 (p + 1)n e + 3(p -1)n f 2 ≤ p < 3 (p + 1)n e + (p -1)(p + 1)n f + 2(p -1)(p -2)n t 3 ≤ p < 4 (p + 1)n e + (p -1)(p + 1)n f + (p -1)(p -2)(p + 1)/2n t 4 ≤ p and Q h := span{(q j ) 1≤j≤n ∂ v }.
The parameter p ∈ N designates the polynomial order of the shape functions. Let us remark that V h,0 represents the usual Raviart-Thomas div-conforming space [START_REF] Raviart | A mixed finite element method for 2-nd order elliptic problems[END_REF]. We define the matrix of the tangential trace operator, denoted by G, as follows

G ij = 1 if w j|Γ = z p i 0 else , with 1 ≤ i ≤ N p z and 1 ≤ j ≤ N p w .
We also define by {•} the vector of degrees of freedom. We will see in Section 8 that a choice on the polynomial order of the unknowns of the system must be made for stability reasons. To this end, let us introduce two parameters (α, β) ∈ N 2 to designate the finite element approximation order.

If T + = Λ 0 -,k -,Z -+ Λ 0,inc -,k -,Z -(the second choice T + = -Λ 0 +,k + ,Z + -Λ 0,inc +,k + ,Z + would only slightly change the formulations), then the discretization of the weak formulation (5.2) is defined in the following proposition. Proposition 6.1. Let us assume that T + = Λ 0 -,k -,Z -+ Λ 0,inc -,k -,Z -. Then, the discretization of the weak formulation (5.2) is given as follows: find (g +,h , g -,h ) ∈ V h,0 × V h,β such that for any test-functions (e 0 i ) i∈{1,...,N 0 z } ∈ V h,0 and (e β j ) j∈{1,...,

N β z } ∈ V h,β ,            Γ h g -,h • e β j dΓ h - Γ h g +,h • e β j dΓ h - Γ h ( a h Z + + m +,h × n Z - ) • e β j dΓ h = 0 Γ h g +,h • e 0 i dΓ h - Γ h g -,h • e 0 i dΓ h + Γ h ( c h Z + + m -,h × n Z - ) • e 0 i dΓ h = 0 , (6.1)
where m +,h is the total exterior magnetic current solution of formulation (6.2) and

m -,h = N α z =1 (G{E -,h }) z α × n,
with E -,h solution of the weak formulation (6.3).

If T -= -Λ sq,Np,θp +,k + ,Z + -Λ sq,inc +,k + ,Z + , then a h and c h are defined through the formulations:

• Find (a h , (Ψ ,h ) ∈{1;...;Np} , (φ ,h ) ∈{1;...;Np} ) ∈ Z h,β × [Z h,β ] Np × [Q h ] Np for any test-functions (z β i ) i∈{1,...,N β z } ∈ Z h,β and (q m ) m∈{1,...,n ∂ v } ∈ Q h ,                                    Γ h C 0 a h • z β i dΓ h - Γ h m +,h × n • z β i dΓ h + Γ h 1 k 2 curl Γ h (m +,h × n) • curl Γ h (z β i ) dΓ h + Np =1 A Γ h ∇ Γ h φ ,h • z β i dΓ h - Γ h 1 k 2 curl Γ h Ψ ,h • curl Γ h (z β i )dΓ h = 0 Γ h Ψ ,h • z β i dΓ h + B Γ h ∇ Γ h φ • z β i dΓ h - Γ h 1 k 2 curl Γ h Ψ ,h • curl Γ h (z β i )dΓ h - Γ h a h • z β i dΓ h = 0 ∀ ∈ {1; ...; N p } Γ h φ • q m dΓ h + Γ h 1 k 2 Ψ • ∇ Γ h q m dΓ h = 0 ∀ ∈ {1; ...; N p } . • Find (c h , (χ ,h ) ∈{1;...;Np} , (ψ ,h ) ∈{1;...;Np} ) ∈ Z h,β × [Z h,β ] Np × [Q h ] Np for any test-functions (z β i ) i∈{1,...,N β z } ∈ Z h,β and (q m ) m∈{1,...,n ∂ v } ∈ Q h ,                                  Γ h C 0 c h • z β i dΓ h - Γ h m -,h × n • z β i dΓ h + Γ h 1 k 2 curl Γ h (m -,h × n) • curl Γ h (z β i ) dΓ h + Np =1 A Γ h ∇ Γ h ψ ,h • z β i dΓ h - Γ h 1 k 2 curl Γ h χ ,h • curl Γ h (z β i )dΓ h = 0 Γ h χ ,h • z β i dΓ h + B Γ h ∇ Γ h ψ • z β i dΓ h - Γ h 1 k 2 curl Γ h χ ,h • curl Γ h (z β i )dΓ h - Γ h c h • z β i dΓ h = 0 ∀ ∈ {1; ...; N p } Γ h ψ • q m dΓ h + Γ h 1 k 2 χ • ∇ Γ h q m dΓ h = 0 ∀ ∈ {1; ...; N p } . If T -= -Λ 0 +,k + ,Z + -Λ 0,inc +,k + ,Z + , then we more simply have: a h = m +,h ×n and c h = m -,h ×n .
The discrete weak formulation of the integral equations (5.3) is defined in the following proposition.

Proposition 6.2. Find (m +,h , j +,h ) ∈ V h,0 × V h,0 such that for any test-function (e 0 i ) i∈{1,...,N 0 z } ∈ V h,0 ,                        1 2 Γ h j +,h • e 0 i dΓ h - Γ h (K k + m +,h Z - ) × n • e 0 i dΓ h + ι Γ h k + Z + Z - (T k + j +,h ) × n •e 0 i dΓ h = 1 2 Γ h g +,h • e 0 i dΓ h + Γ h Z -1 -γ + t E inc × n • e 0 i dΓ h 1 2 Γ h m +,h Z - • e 0 i dΓ h + Γ h (K k + j +,h ) × n • e 0 i dΓ h + ι Γ h k + Z - Z + (T k + m +,h Z - ) × n •e 0 i dΓ h = 1 2 Γ h g +,h × n • e 0 i dΓ h - Γ h γ + t H inc × n • e 0 i dΓ h . (6.2)
Finally, we discretize the variational formulation (5.4) as follows.

Proposition 6.3. Find E -,h ∈ W h,α such that for any test-function

(w i ) i∈{1,...,Nw} ∈ W h,α , Ω h curl E -,h • curl w i dΩ h - Ω h k 2 -E -,h • w i dΩ h + ι Γ h k - Z - Z + c h • Gw i dΓ h = ι Γ h k -Z -g -,h • Gw i dΓ h (6.3)
.

The weak coupling algorithm

In order to have a more compact expression in the discrete setting of the Jacobi algorithm (3.7)-(3.9) at iteration n + 1, we introduce the two families of operators R ± , S ± such that:

{E n+1 + } = R + ({E inc }, {g n + }), {E n+1 -} = R -({0}, {g n -}), {g n+1 -} {g n+1 + } = S -({g n + }, {E n+1 + }) S + ({g n -}, {E n+1 -}) , where {E n+1 -}, {(E n+1 + } and {g n+1 -} {g n+1 + }
are the solutions of the finite element discretizations (6.3), (6.2) and (6.1), respectively. By linearity of the equations, we decompose the field {E n+1 ± } as

{E n+1 ± } = { E n+1 ± } + { En+1 ± },
where

{ E n+1 -} = R -({0}, {0}), { E n+1 + } = R + ({E inc }, {0}), { En+1 ± } = R ± ({0}, {g n ± }).
The quantities { E n+1 ± } are independent of the iteration number n and can hence be written as

{ E ± } := { E n+1 ± }, ∀n ∈ N. Therefore, {g n+1 ± } is such that {g n+1 ± } = S ± ({g n ∓ }, { En+1 ∓ }) + S ± ({0}, { E ∓ }).
One iteration of the weak coupling algorithm then corresponds to

{g n+1 -} {g n+1 + } = A {g n -} {g n + } + b - b + ,
where the quantities

A {g n -} {g n + } are given by A {g n -} {g n + } = S -({g n + }, { En+1 + }) S + ({g n -}, { En+1 -})
.

The information about the incident plane wave is contained in the right-hand side:

b - b + = S -({0}, { E + }) S + ({0}, { E -}. .
The whole weak coupling algorithm can then be recast into a linear system:

(Id -A) {g -} {g + } = b - b +
which can be solved for example by a GMRES iterative solver. Algorithm 1 summarizes the above procedure.

Numerical results

We present in this section some numerical examples to validate and analyze the behaviour of the weak FEM-BEM coupling algorithm 1. All the numerical tests were performed using GmshFEM [START_REF] Royer | Gmsh-Fem: an efficient finite element library based on Gmsh[END_REF], a newly developed open-source finite element library based on Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF] and the open-source boundary element library Bempp-cl [START_REF] Scroggs | Software frameworks for integral equations in electromagnetic scattering based on Calderón identities[END_REF]. Three scatterers are considered: the unit sphere, the unit cube and the same cube with a reentrant corner (see Fig. 2). The objects are illuminated by an electromagnetic plane wave defined by 2. Solve the linear system

E inc (x) = pe ιk + σ inc •x ,
(Id -A) {g -} {g + } = b - b +
by using a Krylov subspace iterative solver (GMRES) [START_REF] Saad | Iterative Methods for Sparse Linear systems[END_REF]. At each iteration, solve the subproblems associated with the systems (6.3) and (6.2).

At convergence, compute the solution

{E ± } = R ± ({E inc }), {0}) + R ± ({0}), {g ± }) = R ± ({E inc }), {g ± }).
where p = (1, 0, 0) T is the polarization vector and σ inc = (0, 0, 1) T is the incidence vector. In the following, we also evaluate the far-field solution through the bistatic Radar Cross Section (RCS) defined by

RCS(σ, σ inc ) := 10log 10 (4π lim r→+∞ r 2 E sc (r, θ, φ) 2 E inc 2 ) (dB),
where σ(θ, φ) = (sin(θ)cos(φ), sin(θ)sin(φ), cos(θ)) T and (r, θ, φ) are the spherical coordinates. For all the presented examples, the GMRES (without restart) iterations are stopped when the initial residual has decreased by a factor of 10 -4 to get an accurate RCS, and we fix the maximum number of iterations to 200. We designate by ( T -; T + ) β α the linear part of the transmission operators selected when considering the finite element approximation order (α; β) as described in Section 6. We define the mesh size h = λ/n λ , with λ := min(λ -; λ + ), setting λ ± = 2π/k ± and where n λ is the number of points per wavelength λ. We always assume that the exterior domain is made up of air so that Z + = Z 0 .

Homogeneous scatterers

Let us start by considering the unit sphere for the parameters: Z -= 1 2 Z + , k -= 2k + and n λ = 5. We report in Fig. 3 the number of GMRES iterations for the weak FEM-BEM coupling vs the exterior wavenumber k + for the operators (-Λ sq,Np,θp +,k + ,Z + ; Λ 0 -,k -,Z -) 3 2 (some comparisons with the choice T + = -Λ 0 +,k + ,Z + will be analyzed later). Several pairs of Padé parameters (N p ; θ p ) are considered. We observe that θ p = π 2 leads to the best convergence and N p = 2 terms are sufficient for the Padé expansion. In the following examples, we will always consider θ p = π 2 and N p = 4. Now, we fix the parameters Z -= 1 2 Z + , k -= 2 and k + = 1, for the unit cube. Fig. 4 shows the influence of the FEM approximation orders (α; β) on the convergence of the weak FEM-BEM coupling with (-Λ sq,Np,θp

+,k + ,Z + ; Λ 0 -,k -,Z -) β α .
The first graph presents the evolution of the number of GMRES iterations in terms of discretization density n λ , where α = 2 and β := 0, ..., 3. We see that for β = 0 and β = 1 the number of iterations strongly increases with n λ . For β = 2 and β = 3 on the other hand, we observe a convergence rate that is independent of n λ , which is an important property for practical problems when refining the mesh. Since the GMRES convergence strongly depends on the spectral properties of the linear part of the weak coupling operator (Id -S π ), the choice of the discretization parameters (α; β) impacts the eigenvalue distribution. To illustrate this point, we report in the second graph the eigenvalues distribution of the linear part of (Id -S π ) in the complex plane when n λ = 10. For the pair (α = 2; β = 0), the convergence is achieved in 62 iterations, which can be explained by the fact that there are eigenvalues with a negative real part. On the other hand, the pair (α = 2; β = 3) leads to an excellent convergence in 9 iterations with a large number of eigenvalues clustering around (1, 0). Since Λ 0 -,k -,Zis a rough approximation of Λ -,k -,Z -, there are also a few eigenvalues distant from (1, 0) but the linear part of (Id -S π ) remains quasi-nilpotent since Λ sq,Np,θp +,k + ,Z + is an accurate approximation of Λ +,k + ,Z + . The numerical results from Fig. 4 suggest that β must fulfill: β ≥ α.

Our goal is now to compare the different transmission operators involved in the formulations, most particularly according to the contrast parameter δ = k -/k + . From now on, we always indicate on the figures the values of the physical parameters k ± and Z ± , and the mesh size n λ . In addition, to simplify the presentation, we introduce the shorter notations: Λ 0 ± := Λ 0 ±,k ± ,Z ± , and Λ sq + := Λ sq,Np,θp +,k + ,Z + , with N p = 4 and θ p = π/2. We present in Fig. 5 the behavior of the weak FEM-BEM coupling for different values of δ. In the case of high-contrast scattering problems, we can clearly see that, as mentioned in Section 3, using T + = -Λ 0 + is more efficient than T + = Λ 0 -. This is even more pronounced when increasing β. To illustrate this property, we set T -= -Λ sq + , with α = 2 and β = 4. We consider an impedance parameter Z x which varies between 10 and 5 × 10 4 . We report in Fig. 6 the evolution of the number of GMRES iterations for the operators

T + = Z -1
x (Id × n) corresponding to various values Z x , for a contrast δ = 13. We observe that the number of iterations is minimal in an interval where Z + seems to be a good approximation of the left point of this interval. Conversely, the value Z -is not a suitable choice since it is out of this range of optimal values. Fig. 7 shows the evolution of the number of GMRES iterations with respect to the exterior wavenumber k + . We observe that (-Λ sq + ; -Λ 0 + ) 3 2 leads to the lowest dependency of the GMRES iterations with respect to the exterior wavenumber. Let us remark that the convergence rates for both (-Λ sq + ; -Λ 0 + ) 3 2 and (-Λ sq + ; Λ 0 -) 3 2 are quasi-identical since the contrast material δ = 2 is low. As it can be seen from Figures 5 and7, the increase of the number of iterations is more significant for the cube and the cube with reentrant corner than for the sphere. This is due to the fact that the approximations (4.1) are a priori formally derived for smooth and convex objects. Fig. 8 shows the evolution of the number of GMRES iterations according to the discretization density n λ . For (-Λ sq + ; ∓Λ 0 ± ) 3 2 , the rate of convergence is mesh independent while it clearly deteriorates with (-Λ 0 + ; ∓Λ 0 ± ) 3 2 . From these results, we conclude that the weak FEM-BEM coupling with the choice of operators (-Λ sq + ; -Λ 0 + ) β α , for β > α, is very robust and well-suited for solving high-frequency and high-contrast scattering problems.

In order to validate the weak FEM-BEM coupling, the relative errors in L 2 t (Γ)-norm between the analytical (via a Mie series expansion [START_REF] Born | Principle of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light[END_REF]) and numerical electric/magnetic currents are computed in the sphere case for different densities n λ (see Fig. 9). We designate by e γ - t H and e γ - t E the relative errors associated with the electric and magnetic currents, respectively. We remark a lower accuracy of the magnetic current for the cube and the cube with reentrant corner. However, this loss of accuracy can easily be solved by considering the interior problem (2.1a) with known numerical surface electric current during a post-processing step. In Fig. 10, we draw RCS(σ(θ, 0), σ inc ) as a function of the angle θ (i.e. in the plane {x 2 = 0}) for different GMRES tolerances and for the recommended operators (-Λ sq + ; -Λ 0 + ) 3 2 . For the sphere, the reference analytical far-field [START_REF] Born | Principle of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light[END_REF] is also computed. As expected, the results confirm that the RCS converges when the GMRES tolerance decreases and that a tolerance of 10 -4 is enough to get an accurate RCS.

(-Λ 0 + ; Λ 0 -) 3 2 (-Λ sq + ; Λ 0 -) 3 2 (-Λ sq + ; Λ 0 -) 4 2 (-Λ 0 + ; -Λ 0 + ) 3 2 (-Λ sq + ; -Λ 0 + ) 3 2 (-Λ sq + ; -Λ 0 + ) 4 2 0 5 10 15 20
10 1 10 2 δ Number of GMRES iterations Cube (-Λ 0 + ; Λ 0 -) 3 2 (-Λ sq + ; Λ 0 -) 3 2 (-Λ sq + ; Λ 0 -) 4 2 (-Λ 0 + ; -Λ 0 + ) 3 2 (-Λ sq + ; -Λ 0 + ) 3 2 (-Λ sq + ; -Λ 0 + ) 4 2 0 5 10 15

Inhomogeneous scatterers

To end, we analyze the case of inhomogeneous scatterers. More specifically, we choose k -and Z -as the space dependent functions given by

k -= δk + e -x 2 , Z -= Z + δ e x 2 ,
with contrast parameter δ ∈ R * + . The values of k + , δ and n λ are indicated on the figures. We present in Fig. 11 the behavior of the weak FEM-BEM coupling for different contrast parameters δ. In Fig. 12, we report the number of GMRES iterations required for solving the weak FEM-BEM coupling as a function of the exterior wavenumber k + . Fig. 13 shows the effect of the mesh refinement n λ over the number of GMRES iterations. Finally, we plot on Fig. 14 the bistatic RCS obtained with various GMRES tolerances for the operators (-Λ sq + ; -Λ 0 + ) 3 2 . As can be seen on these different figures, we conclude that we obtain results similar to the ones for the homogeneous case.

Conclusion

In this work, we proposed a new well-conditioned weak coupling of boundary element and highorder finite element methods for time-harmonic electromagnetic scattering. The weak coupling was formulated as a non-overlapping Schwarz domain decomposition method, where the transmission conditions are constructed through Padé localized approximations of the Magnetic-to-Electric operators. Numerical results validating the new approach are presented for both homogeneous and inhomogeneous obstacles. The number of iterations required to solve the weak FEM-BEM coupling is only slightly dependent on the frequency, the mesh size and the contrast coefficient. Extensions to scatterers with partially coated boundaries and coupling with H-matrix solvers are currently underway.

Figure 1 :

 1 Figure 1: The electromagnetic transmission-scattering problem.

Proposition 5 . 5 . 1 2

 551 Denoting by ; the duality pairing u; v = Γ u×n•v dΓ between H -

Algorithm 1 :

 1 Weak coupling algorithm with Krylov solver. 1. Compute the right-hand side b - b + .

Figure 2 :

 2 Figure 2: Scatterer geometries: sphere, cube and cube with reentrant corner.

6 Figure 3 :

 63 Figure 3: Influence of parameters Np and θp on the convergence of the weak FEM-BEM coupling.

Figure 4 :

 4 Figure 4: Influence of parameters (α;β) on the convergence of the weak FEM-BEM coupling.

2 δ

 2 Number of GMRES iterationsSphere

2 Figure 5 :Figure 6 :

 256 Figure 5: Number of GMRES iterations for the weak FEM-BEM coupling vs the contrast parameter δ (with k-= δk+, Z-= Z+/δ, n λ = 5).

2 Figure 7 :Figure 8 :

 278 Figure 7: Number of GMRES iterations for the weak FEM-BEM coupling vs the exterior wavenumber k+ (with k-= 2k+, Z-= Z+/2, n λ = 5).

2 Figure 9 :

 29 Figure 9: Relative errors between the analytical and numerical magnetic/electric surface currents as functions of the discretization density n λ (with k-= 4, k+ = 2, Z-= Z+/2).

- 1 tol = 10 - 2 tol = 10 - 3 tol = 10 - 4 Figure 10 : 2 Figure 11 : 2 Figure 13 : 1 tol = 10 - 2 tol = 10 - 3 tol = 10 - 4 Figure 14 :

 110210310410211213110210310414 Figure 10: Convergence of the RCS with respect to the GMRES tolerance for the weak FEM-BEM coupling (with k-= 6, k-= 3, Z-= Z+/2, n λ = 10).

  where the fields (E -; H -) ∈ H(curl , Ω -) 2 and (E + ; H + ) ∈ H loc (curl , Ω + ) 2 are the solutions of the respective boundary-value problems

t (div Γ , Γ), and ∀x ∈ R

\ Γ, T k j(x) = Γ G k (x -y)j(y) dy + 1 k 2 ∇ Γ G k (x -y)div Γ j(y) dy

Acknowledgments. This work was supported by Thales Defence Mission Systems (DMS) through the CIFRE contracts 2015/0197 and 2018/1609. This work was funded in part through an ARC grant for Concerted Research Actions (ARC WAVES 15/19-03), financed by the Wallonia-Brussels Federation of Belgium. Computational resources have been provided by the Consortium des


Eugene M Burreson

Michel Segonzac
email: michel.segonzac@ifremer.fr.

Morphological variability of Bathybdella sawyeri (Hirudinida: Piscicolidae) from hydrothermal vents on the Galápagos Rift and the South East Pacific Rise

Keywords: Hirudinida, Piscicolidae, East Pacific Rise, hydrothermal vent, Bathybdella, Bythograeidae, Cyanagraea

The deep-sea fish leech Bathybdella sawyeri Burreson has recently been collected from two sites on the South East Pacific Rise at approximately14° and 17°S, 2623-2586 m depth. Leeches were obtained from washings of vent mussels Bathymodiolus thermophilus (Mytilidae) and clams Calyptogena magnifica (Vesicomyidae) collected with a grab by a submersible, or were found attached to the cheliped of the crab Cyanagraea praedator (Bythograeidae) collected in baited traps at 14°S. Most individual leeches had morphology typical of B. sawyeri as originally described, with suckers of approximately equal size, but all specimens collected from the crab C. praedator had unusually large oral suckers that were distinctly larger than the caudal sucker. However, internal morphology of the specimens from C. praedator was identical to B. sawyeri on the basis of serial sections; thus all specimens were identified as B. sawyeri. The only previous report of B. sawyeri is from the Galápagos Rift, so these collections extend the range to the South East Pacific Rise.

Introduction

Bathybdella sawyeri Burreson was described from material collected on the first Galápagos Rift Biology Expedition in the eastern Pacific Ocean in 1979 [START_REF] Burreson | A new deep-sea leech, Bathybdella sawyeri, n. gen., n. sp., from thermal vent areas on the Galápagos Rift[END_REF]. Most specimens collected have been found free in washings from clumps of mussels Bathymodiolus thermophilus Kenk & Wilson (Mytiloida: Mytilidae), the clam Calyptogena magnifica Boss & Terner (Veneroida: Vesicomyidae), or the vestimentiferan tube worm Riftia pachyptila Jones (Sabellida: Siboglinidae). Bathybdella sawyeri is unusual in a number of respects including depth of capture, high abundance, and reproductive system anatomy. Like all members of the family Piscicolidae, the leech feeds on fish blood as evidenced by nucleated red blood cells in the gut, although no specimen has ever been collected from a fish host. Because of its common association with invertebrates, it seems likely that B. sawyeri leaves the fish host after each blood meal and attaches to the hard surface of various invertebrates or other hard substrate until ready to feed again.

Recent collections at two sites along the South East Pacific Rise have extended the range of B. sawyeri and also provided specimens from the crab Cyanagraea praedator de Saint Laurent (Decapoda; Bythograeidae) that differ from B. sawyeri in size of the oral sucker. In this paper we describe the morphological variability and discuss some aspects of the biology of the leech.

Material and Methods

The specimens were collected during the French diving cruise BIOSPEEDO, using R/V L'Atalante and D/S Nautile, from 31 March to 13 May 2004, Chief scientist D. Jollivet (CNRS, Roscoff), and carried out along the South East Pacific Rise between 7°24'S and 21°33'S [START_REF] Jollivet | The BIOSPEEDO Cruise; a new survey of hydrothermal vents along the South East Pacific Rise from 7°24'S to 21°33[END_REF]) (Fig. 1). The majority of the leech specimens were found free from washings of the bivalves C. magnifica and B. thermophilus or the crabs Bythograea thermydron Williams (Decapoda: Bythograeidae) and C. praedator and of what is likely a new genus and species of fish (Ophidiiformes: Bythitidae) collected using the grab of the submarine or with baited traps. Eleven specimens were attached to the cheliped of the crab C. praedator collected in baited traps, and two were on the external part of the hinge of a clam C. magnifica. Seven specimens came from four dives on the Oasis site (17°25.38'S, 113°12.29'W, 2586 m) and eleven specimens from the Pagoda site (13°59.06'S, 112°28.94'W, 2623 m), nearly 6 km N of Oasis. Onboard, the leeches were fixed either in formalin or in 95% ethanol for DNA analysis.

Material examined

Oasis site, 17°25.38'S, 113°12.29'W, 258617°25.38'S, 113°12.29'W, m. PL 158217°25.38'S, 113°12.29'W, , 22.04.2004, baited , baited In addition to the newly collected specimens, the holotype (USNM 65773) and paratypes (USNM 65774, 65775) of B. sawyeri housed in the U. S. National Museum were re-examined, as well as 109 individuals collected on Woods Hole Oceanographic Institution Alvin Dive 884. That dive occurred on January 25, 1979 on the Galápagos Rift (00°47.69'N, 86°07.7'W). Leeches were from washings of tubeworms R. pachyptila and bivalves B. thermophilus. These leeches are housed at the U. S. National Museum, but are uncatalogued.

Habitat of the hosts

Oasis vent site

The Oasis site is characterized by a milky fluid diffusing from crevices, where the temperature varies between 2 and 7°C, surrounded by a very abundant vent community composed of mussels, B. thermophilus, clams, C. magnifica, a stalked barnacle, Neolepas n. sp., sea anemone Chondrophellia cf. coronata (Verrill) (Actinostolidae), the crabs B. thermydron and Munidopsis subsquamosa Henderson (Galatheidae), and a recently described shrimp, Nematocarcinus burukovskyi Komai & Segonzac (Nematocarcinidae) [START_REF] Komai | Two new species of Nematocarcinus A. Milne-Edwards, 1881 (Crustacea: Decapoda: Caridea: Nematocarcinidae) from hydrothermal vents on the North and South East Pacific Rise[END_REF]. Numerous fish occur in this environment such as Ilyophis saldanhai Karmovskaya & Parin (Synaphobranchidae), Thermychthys hollisi (Cohen, Rosenblatt & Moser) (Bythitidae) (Nielsen & Cohen in press), a very abundant Ophidiidae species, and one hagfish species.

Pagoda vent site (BS 17)

This site is composed of two black smokers covered by an apparently declining tubeworm polychaete colony (Alvinellidae); temperature varied from 2°C in the vicinity of the black smokers to 10°C near the tubeworm colony. A surprising abundance of the crab C. praedator (12 specimens were collected in a baited trap) was observed among the alvinellid colony. A small number of fish (Ophidiidae) was present within 3 m of the trap.

Results

All leeches from the Oasis site were typical B. sawyeri with maximum length of 11 mm and suckers of equal size (maximum 1.0 mm in diameter). Photographs taken in life showed them to be colorless except for the red blood in the gut that was visible through the transparent body wall. Most leeches from Oasis were engorged with blood (Fig. 2).

All leeches from Pagoda site were attached to the crab C. praedator. They were unusual in that the overall length of the leeches was larger (13 mm) than typical B. sawyeri and the oral sucker was distinctly larger than the caudal sucker (Fig. 3). There was little variation in the leeches from Pagoda; they were all 13 mm in length and all had an oral sucker with a diameter of 1.3 mm and a caudal sucker with a diameter of 1.0 mm. A single leech from Pagoda that was serial sectioned transversely for its entire length revealed an internal anatomy identical with B. sawyeri in all respects.

All 109 leeches examined from Alvin dive 884 were typical B. sawyeri in terms of sucker size, with one possible exception. A single leech, 9 mm in length, had an oral sucker that was larger than the caudal sucker, 0.8 mm vs 0.6 mm. The other leeches were remarkably consistent with the holotype and paratypes of B. sawyeri and had suckers that were approximately equal in size.

Discussion

The collections reported here extend the range of B. sawyeri from the Galápagos Rift to the South East Pacific Rise, and suggest that the leech is widely distributed at vent sites in the Pacific Ocean. None of the leeches from the Oasis site was sectioned, but externally they are identical to B. sawyeri. The leeches from the Pagoda site were not originally identified as B. sawyeri because of the large oral sucker. However, the serial sections revealed that the internal anatomy of these leeches is identical to B. sawyeri, including 6 pairs of testisacs, and the very unusual reproductive system that has a large bilobed sprematheca connected via conducting tissue to paired, deep lateral invaginations of the epidermis [START_REF] Burreson | A new deep-sea leech, Bathybdella sawyeri, n. gen., n. sp., from thermal vent areas on the Galápagos Rift[END_REF]. These leeches clearly belong in the genus Bathybdella based on internal anatomy, but the oral sucker is unusually large for B. sawyeri. It is possible that the leeches from the Pagoda site are a separate species of Bathybdella; however, we have chosen to identify the leeches with the large oral sucker as B. sawyeri because the internal anatomy is identical with that species. DNA sequence data will likely be necessary to confirm the identification or provide evidence that it is a separate species.

It is also interesting that all the leeches from Pagoda were collected from the crab C. praedator, and none of the leeches from Oasis were collected on this crab. The association with C. praedator is likely not specific, however. Many marine leeches are known to associate with decapod crustaceans, mainly for hard surface attachment sites after they leave the fish host, and often for sites of cocoon deposition [START_REF] Meyer | Leeches symbiotic on Arthropoda, especially decapod crustacea[END_REF][START_REF] Daniels | The biology of the leech Myzobdella lugubris infesting blue crabs and catfish[END_REF][START_REF] Sloan | Cocoon deposition on three crab species and fish parasitism by the leech Notostomum cyclostoma from deep fjords in northern British Columbia[END_REF][START_REF] Khan | Life cycle studies on Arcto-boreal leeches (Hirudinea)[END_REF].

Leeches in the collections reported here were also recovered from clams, supporting the notion that they are attaching to any hard substrate available after leaving the fish host.

Specimens of B. sawyeri never have been recovered from a fish host, but few fish have actually been collected from vent sites. Thus, it is tempting to speculate that the leech feeds on the hemolymph of the clam C. magnifica or the tube worm R. pachyptila. Both these species have hemoglobin free in the hemolymph giving it a red color resembling the red color of the leech gut contents. However, histological sections of the leech gut reveal the abundance of nucleated red blood cells typical in size and morphology to fish red blood cells. Their presence is evidence that B. sawyeri, like all other leeches in the family Piscicolidae, feeds on fish blood. Although they were not directly observed on the two Bythitidae fish collected in the baited trap C2 from the Oasis site, these leeches are most probably parasitic of this fish. They may, however, be parasitic also on other fish species present at the sites visited during the BIOSPEEDO cruise between 7°24'S and 21°33'S [START_REF] Jollivet | The BIOSPEEDO Cruise; a new survey of hydrothermal vents along the South East Pacific Rise from 7°24'S to 21°33[END_REF]) including a species of the family Ophidiidae, visible in large numbers on the video documents taken at the same sites, but not captured. Marine leeches vary considerably in their host specificity. Many species occur only on a single species of fish [START_REF] Burreson | Two new species of Malmiana (Hirudinea: Piscicolinae) from Oregon coastal waters[END_REF], while other species feed on a wide variety of fish hosts [START_REF] Burreson | The role of host biology, vector biology, and temperature in the distribution of Trypanoplasma bullocki infections in the lower Chesapeake Bay[END_REF]; the host specificity of B. sawyeri is unknown.

The unusually high abundance of B. sawyeri at vent sites, noted previously by [START_REF] Burreson | A new deep-sea leech, Bathybdella sawyeri, n. gen., n. sp., from thermal vent areas on the Galápagos Rift[END_REF], is supported by the collection of over 100 specimens from a single Alvin dive at the Galápagos Rift in 1979. The high abundance of leeches observed on or among invertebrates suggests that B. sawyeri leaves the fish host after each blood meal and attaches to any hard substrate available. When ready to feed again individuals must reattach to a fish, but finding a fish host (by random contact) is likely not difficult considering the high abundance of demersal fishes in the vent areas. When ready to reproduce, marine leeches mate and deposit cocoons (egg cases), containing usually a single egg, on any hard substrate available, including the carapaces of decapod crustaceans. A search for clusters of leech cocoons on the hard surface of the various invertebrates or other hard substrate at thermal vent sites should be productive. Almost nothing is known of the biology of B. sawyeri at vent sites, but the high abundance of the leech suggests that the life cycle of B. sawyeri is completed in the vicinity of the vents, with fish hosts for feeding, and hard substrate for cocoon deposition, in abundance.

FIGURE 1 .

 1 FIGURE 1. Map of collection sites (open circles) on the Galápagos Rift and South East Pacific Rise and the known distribution of B. sawyeri (closed circles).

FIGURES 2- 3 .

 3 FIGURES 2-3. Bathybdella sawyeri. 2: specimens from Oasis showing lack of pigmentation and the gut full of blood, ventral view. 3: specimen from Pagoda taken from the crab Cyanagraea praedator, lateral view. o, oral sucker. Scale bars = 2.0 mm.

  MNHN AH1-13VV. PL 1586, 27.04.2004, basket 1: 2 sps. in formalin, with mytilid B. thermophilus, MNHN AH1-14VV. PL 1589, 30.04.2004, baited trap B3: 3 sps. in formalin, with 5 crabs B. thermydron, MNHN AH1-15VV. Pagoda site (marker BS 17), 13°59.06'S, 112°28.94'W, 2623 m. PL 1592, 04.05.2004, baited trap F2: 6 sps. (2 in formalin, 4 in ethanol), on crab C. praedator, MNHN AH1-16VV. Baited trap F2: 1 sp. in formalin, serial transverse sections, on crab C. praedator. Basket 1: 4 sps. in formalin, on crab C. praedator, MNHN AH1-17VV.

	trap C2: 1 specimen in formalin, found with 3 crabs B. thermydron and 2
	bythitid fish gen. et sp. nov. Muséum National d'Histoire Naturelle (MNHN), Paris, France AH1-12VV.
	PL 1583, 24.04.2004, basket 1: 1 sp. in formalin, with clams C. magnifica,

Acknowledgements

We thank Didier Jollivet (Chief scientist of the BIOSPEEDO cruise, Roscoff, France) for collecting the specimens studied here. Also thanks to Violaine Martin (Ifremer) for Fig. 1 and Stéphane Hourdez (CNRS, Roscoff) for the pictures of B. sawyeri. The cruise was supported by the French programme DORSALES and the Ecchis GDR (CNRS and Ifremer).