Marc Bonnet
email: bonnet@lms.polytechnique.fr

Inverse acoustic scattering by small-obstacle expansion of misfit function

This article concerns an extension of the topological derivative concept for 3D inverse acoustic scattering problems, whereby the featured cost function J is expanded in powers of the characteristic size ε of a sound-hard scatterer about ε = 0. The O(ε 6 ) approximation of J is established for a small scatterer of arbitrary shape of given location embedded in an arbitrary acoustic domain, and generalized to several such scatterers. Simpler and more explicit versions of this result are obtained for a centrally-symmetric scatterer and a spherical scatterer. An approximate and computationally fast global search procedure is proposed, where the location and size of the unknown scatterer is estimated by minimizing the O(ε 6 ) approximation of J over a search grid. Its usefulness is demonstrated on numerical experiments, where the identification of a spherical, ellipsoidal or banana-shaped scatterer embedded in a acoustic half-space from known acoustic pressure on the surface is considered.

Introduction

The reconstruction of obstacles embedded in a three-dimensional acoustic medium [START_REF] Colton | Regularized sampling method for solving threedimensional inverse scattering problems[END_REF][START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF] is a challenging subject with applications to e.g. sonar detection and medical imaging, and also provides useful insight into e.g. elastodynamic or electromagnetic inverse scattering. This class of non-linear inverse problem may be formulated in terms of the minimization of a cost function featuring the experimental data and (possibly) prior information. Such cost functions are non-convex and exhibit local minima. Global search techniques, e.g. evolutionary algorithms [START_REF] Michalewicz | How to solve it: modern heuristics[END_REF] or parameter-space sampling methods [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF] require large numbers of cost functions evaluations, which makes them extremely expensive in the present context due to the high computational cost entailed by each forward scattering solution. Hence, traditional minimization methods [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF][START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized Newton method[END_REF][START_REF] Nishimura | Determination of cracks having arbitrary shape with boundary integral method[END_REF], or Newton-type algorithms for solving the observation equations [START_REF] Kress | A hybrid method for two-dimensional crack reconstruction[END_REF][START_REF] Potthast | On the convergence of a new Newton-type method in inverse scattering[END_REF], are often used, as they may converge within a moderate number of forward solutions if the obstacle is described in terms of a small number of parameters. To optimize computational efficiency, such solution techniques are used in conjunction with shape sensitivity techniques [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF][START_REF] Kirsch | The domain derivative and two applications in inverse scattering theory[END_REF][START_REF] Nintcheu Fata | A computational basis for elastodynamic cavity identification in a semi-infinite solid[END_REF][START_REF] Sokolowski | Introduction to shape optimization[END_REF] or level-set methods [START_REF] Litman | Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set[END_REF].

Still, the stand-alone use of gradient-based algorithms for such purposes is not always satisfactory due to their strong dependence on reliable prior information about the geometry of the hidden object(s). This has prompted the development of alternative, non-iterative, approaches, which may be used either in isolation or as a preliminary step providing adequate initial guesses to subsequent iterative schemes. Now referred to as "sampling" or "probe" methods [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF], they include in particular the linear sampling method [7, 9-11, 24, 32]. Another way of sampling spatial regions that may contain unknown objects, to which this article contributes, exploits the concept of topological sensitivity, whereby the perturbation of a cost function with respect to the creation of a scatterer B ε (a) of given shape, location a, and small characteristic radius ε, is quantified as a function of a. Letting J(ε; a) denote the value achieved by the cost function used for solving the inverse problem for the trial scatterer B ε (a), then in 3D situations the topological derivative T 3 (a) appears through an expansion of the form

J(ε; a) = J(0) + ε 3 T 3 (a) + o(ε 3 ) (1)
The concept of topological derivative was first used for the topological optimization of mechanical structures [START_REF] Eschenauer | Bubble method for topology and shape optimization of structures[END_REF][START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Lochpositionierungskriterien[END_REF] by means of algorithms where "excess" material is iteratively removed until a satisfactory shape and topology is reached [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF]. More recently, other investigations have studied the topological sensitivity as a preliminary sampling tool for inverse scattering problems, providing estimates of location, size and number of defects which can then (for example) be used as initial guesses in subsequent minimization-based inversion procedures. Previous contributions in this direction address wave-based identification of cavities or penetrable inclusions in 3D elastic solids using frequency-domain [START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF][START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] or time-domain [START_REF] Bonnet | Topological sensitivity for 3D elastodynamic and acoustic inverse scattering in the time domain[END_REF][START_REF] Dominguez | Time domain topological gradient and time reversal analogy: an inverse method for ultrasonic target detection[END_REF] formulations, and of penetrable inclusions in 3D acoustic media [START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF]. Other closely related investigations include [START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF][START_REF] Jackowska-Strumillo | Topological optimization and inverse problems[END_REF] for 2D elastostatics, [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] for 2D linear acoustics, and the comprehensive study of small-inclusion asymptotics of [START_REF] Ammari | Reconstruction of small inhomogeneities from boundary measurements[END_REF].

In this article, building on [START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF][START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], an extension of the topological derivative is proposed whereby J(ε; a) is expanded further in powers of ε. Specifically, the expansion to order O(ε 6 ) for misfit functions associated with 3D acoustic scattering by a soundhard obstacle of size ε embedded in a medium occupying a domain of arbitrary shape is established. The chosen order O(ε 6 ) stems from the fact that, for misfit functions J of least-squares format, the perturbations of the residuals featured in J are of order O(ε 3 ) under the present conditions. The expansion will be found to have the form J(ε; a) = J 6 (ε; a) + o(ε 6 ), J 6 (ε; a) ≡ J(0) + T 3 (a)ε 3 + T 4 (a)ε 4 + T 5 (a)ε 5 + T 6 (a)ε 6 (2) where coefficients T 3 (a), T 4 (a), T 5 (a), T 6 (a) depend on the assumed shape of the scatterer and are expressed in terms of the free field and the adjoint field, i.e. the response of the reference medium to the probing excitation and to the adjoint excitation (defined in terms of the misfit function density), and also (for T 6 (a)) on the leading O(ε 3 ) contribution to the scattered field at measurement location. The functions T 3 (a), . . . , T 6 (a) can be computed for a spanning a 3-D search grid G at a computational cost which is of the order of at most a small number of forward solutions in the reference medium. Minimizing the polynomial J 6 (ε; a) is a simple and inexpensive task, which can therefore be performed for a dense search grid G, thereby defining an approximate, computationally fast, global search procedure. The values of a and ε leading to a global minimum of J 6 (ε; a) over a ∈ G can then be used as either a stand-alone estimate of the sought scatterer or an initial guess for a subsequent refined search procedure. The usefulness of this proposed exploitation of expansion J 6 (ε; a) will be demonstrated through numerical experiments.

This article is organized as follows. The forward and inverse problems of interest are defined in Section 2. Detailed expressions for the coefficients of expansion (2) are established in Section 5 for a scatterer of arbitrary shape embedded in an arbitrary domain, based on a methodology whose main components consist of an adjoint-solution framework (Section 4) and an expansion of the scattered field on the obstacle surface (Section 3). The useful special case of a centrally-symmetric scatterer is further shown in Section 5.2 to lead to significantly simpler formulae, which moreover become explicit for spherical scatterers (Section 5.3). The generalization to several scatterers is treated next (Section 6). After discussing computational issues and links to other approaches in Section 7, a simple approximate global search procedure based on J 6 (ε; a) is proposed and demonstrated on numerical experiments in Section 8.

Forward and inverse problems

Consider a three-dimensional domain Ω, either bounded or unbounded, filled with an acoustic medium characterized by wave velocity c and mass density ρ; this configuration will be referred to as the reference (i.e. obstacle-free) medium. An unknown sound-hard obstacle B true , bounded by the closed surface Γ true , is embedded in the reference medium, so that the acoustic region surrounding the obstacle is Ω true = Ω \ (B true ∪ Γ true ). For identification purposes, the obstacle is illuminated by a free (i.e. incident) field defined by the boundary-value problem

(∆ + k 2 )u = 0 (in Ω), p = p D (on S N ), u = u D (on S D ), (3)
where the boundary S of Ω is divided into complementary subsets S N and S D supporting Neumann and Dirichlet data p D and u D , k = ω/c is the wavenumber, n is the normal on S ∪ Γ true outward to Ω true , and with p ≡ ∇u•n. For simplicity, it is assumed that ω is not an eigenfrequency of any of the boundary-value problems arising in this work. Considering a trial obstacle B bounded by Γ , the prescribed excitation (p D , u D ) gives rise in Ω = Ω \ (B ∪ Γ ) to an acoustic field u which can be conveniently be decomposed as

u = u + v , (4)
where v , the scattered field in the region Ω = Ω \ (B ∪ Γ ) surrounding B , solves

(∆ + k 2 )v = 0 (in Ω ), q = 0 (on S N ), v = 0 (on S D ), q = -p (on Γ ), (5)
with q ≡ ∇v •n.

The inverse problem considered here consists in identifying the unknown scatterer B true from supplementary data consisting of known values u obs of the acoustic pressure and -(1/ρ)p obs of the normal wall acceleration, collected respectively on the measurement surfaces S obs N ⊂ S N and S obs D ⊂ S D . The misfit between observations u obs , p obs and their acoustic predictions u , p ≡ ∇u • n for a trial obstacle B true is quantified through a cost function J (Ω ) to be minimized. Generic cost function having the format

J (Ω ) = S N ϕ N u R (ξ), u I (ξ), ξ dΓ + S D ϕ D p R (ξ), p I (ξ), ξ dΓ (6)
are considered, where ϕ N and ϕ D are C 2 functions with respect to their first two arguments and the subscripts 'R' and 'I' indicate the real and imaginary parts of a complex number (i.e. w R = Re(w), w I = Im(w)). With adequate definitions for ϕ N and ϕ D , the format (6) allows to accommodate partial (or even empty) measurement regions.

For instance, a weighted output least-squares cost function associated to measurement u obs on S obs N ⊂ S N (with non-negative weighting function W (ξ) > 0), commonly used for such purposes, corresponds to

ϕ(w R , w I , ξ) = 1 2 W (ξ) w(ξ)-u obs (ξ) 2 1 S obs N . (7)
In what follows, attention will focus on the case of trial obstacles of small size ε and given location a and shape. The main objective of this article is to establish an expansion of cost functions of format [START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF] with respect to ε, whose coefficients depend on a, and formulate an approximate global search method which exploits that dependence.

Adjoint solution framework for expansion of misfit function

Preliminaries

Let B ⊂ R 3 is a fixed bounded open set containing the origin, with boundary S and volume |B|. With this definition, one may consider the introduction of a small soundhard obstacle B ε (a) = a + εB, with boundary Γ ε (a) = a + εS , of size ε > 0 and centered at a given point a (thereafter referred to as a sampling point). The acoustic region surrounding the small obstacle is then Ω ε (a) = Ω \ (B ε (a) ∪ Γ ε (a)).

To investigate small-obstacle approximations of cost functions (6), the trial domain Ω is now defined in terms of a trial obstacle B = B ε (a) of small size ε, i.e. Ω = Ω ε (a). Accordingly, let u ε (•; a) = u(•)+v ε (•; a) denote the solution to the scattering problem (4), [START_REF] Bonnet | Topological sensitivity for 3D elastodynamic and acoustic inverse scattering in the time domain[END_REF] with Ω = Ω ε (a), and define J(ε; a) by

J(ε; a) = J Ω ε (a) = S N ϕ N u ε R (ξ), u ε I (ξ), ξ dΓ ξ + S D ϕ D p ε R (ξ), p ε I (ξ), ξ dΓ ξ , (8)
with p ε ≡ ∇u ε •n. For notational convenience, explicit references to a will sometimes be omitted in the sequel, e.g. by writing J(ε) or u ε (ξ) instead of J(ε; a) or u ε (ξ; a).

Expansion of misfit function using adjoint solution

Misfit functions used in applications often are of least-squares format, i.e. depend quadratically on the measurement residuals. With this in mind, the desired polynomial approximation of J(ε) is sought by exploiting an expansion of (8) to second order in (v ε ), i.e.:

J(ε) = J(0) + S N Re ϕ N,u v ε dΓ + S D Re ϕ D,p q ε dΓ + 1 2 a,b=R,I S N ϕ N,ab v ε a v ε b dΓ + S D ϕ D,ab q ε a q ε b dΓ + o |v ε | 2 L 2 (S N ) , |q ε | 2 L 2 (S D ) , (9)
(which is exact, i.e. has a zero remainder, for least-squares misfit functions), having set q ε ≡ ∇v ε •n and

ϕ N,u = ∂ϕ N ∂u ε R -i ∂ϕ N ∂u ε I u ε =u , ϕ D,p = ∂ϕ D ∂p ε R -i ∂ϕ D ∂p ε I p ε =p , (10)
ϕ N,ab = ∂ 2 ϕ N ∂u ε a ∂u ε b u ε =u , ϕ D,ab = ∂ 2 ϕ D ∂p ε a ∂p ε b p ε =p (a, b = R, I). (11)
Let the adjoint field û be defined as the solution of the adjoint problem

(∆ + k 2 )û = 0 (in Ω) , p = ϕ N,u (on S N ) , û = -ϕ D,p (on S D ). (12)
(with p ≡ ∇û•n). Then, the third Green's formula applied to the fields û and v ε on the domain Ω ε leads, by virtue of the boundary conditions in [START_REF] Bonnet | Topological sensitivity for 3D elastodynamic and acoustic inverse scattering in the time domain[END_REF] and [START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF], to the identity

S N ϕ N,u v ε dΓ + S D ϕ D,p q ε dΓ = - Γε û p dΓ - Γε p v ε dΓ = 0 (13)
Applying the divergence formula to the third integral of (13) and using the resulting equality in [START_REF] Colton | Regularized sampling method for solving threedimensional inverse scattering problems[END_REF], one arrives at

J(ε) = J(0) + Re Bε ∇u•∇û -k 2 uû dV - Γε p v ε dΓ + 1 2 a,b=R,I S N ϕ N,ab v ε a v ε b dΓ + S D ϕ D,ab q ε a q ε b dΓ + o(|v ε | 2 L 2 (S N ) , |q ε | 2 L 2 (S D ) ). ( 14
)
Summary of previous results on topological sensitivity. The leading contribution to J(ε) can be determined on the basis of identity ( 14) truncated to first order in (v ε , q ε ) (i.e. without the last two integrals). It has been found in previous studies [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF][START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF] to have the form

J(ε; a) = J(0) + ε 3 T 3 (a) + o(ε 3 ) ( 15
)
in terms of the topological derivative T 3 (a), given in the present context of acoustic scattering by sound-hard obstacles by

T 3 (a) = Re ∇û•A 11 •∇u -|B| k 2 ûu (a) (16)
where the second-order 'polarization tensor' A 11 has been established for any obstacle shape S (including the more general case of penetrable inclusions) in e.g. [START_REF] Ammari | Recovery of small inhomogeneities from the scattering amplitude at a fixed frequency[END_REF][START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF][START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF].

For the simplest case of a sound-hard spherical obstacle, where B is the unit sphere, one has the explicit expression

A 11 = 2πI. ( 17
)
(with I the second-order identity tensor). Moreover, the leading asymptotic behaviour of the scattered field is characterized by

v ε (x) = ε 3 W (x) + o(ε 3 ), q ε (x) = ε 3 Q(x) + o(ε 3 ) (x ∈ S) (18)
(having set Q(x) = ∇W (x)•n(x)) on the external surface, and

v ε (x) = εV 1 (x -a)/ε + o(ε) (x ∈ Γ ε ) (19)
on the surface of the small scatterer, where the functions W and V 1 are known and depend on the obstacle shape S (see equations ( 61) and ( 41)).

Derivation of expansion of J(ε): methodology and notation

To include the leading contribution as ε → 0 of the quadratic terms v ε a v ε b and q ε a q ε b , expansion of J(ε) must, in view of ( 14) and [START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized Newton method[END_REF], be performed up to at least order O(ε 6 ). Since [START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF] involves integrals over vanishing supports B ε and Γ ε , this task is facilitated by scaling the position vector ξ in B ε and on Γ ε according to:

ξ = a + ε ξ (ξ ∈ B ε , ξ ∈ B) (20)
In particular, this mapping recasts integrals over B ε and Γ ε into integrals over B and S , respectively, and transforms the differential volume and area elements according to

(a) dV ξ = ε 3 d Vξ (ξ ∈ B ε , ξ ∈ B), (b) dΓ ξ = ε 2 d Γξ (ξ ∈ Γ ε , ξ ∈ S ) (21)
Without loss of generality, a can be chosen as the center of B ε , i.e. such that

B ξ d Vξ = 0. (22)
The following simple integral identities are collected for later reference: where (23a) and (23c) exploit [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] and the minus sign in (23b) stems from the fact that n is the inward unit normal to S .

In view of (21b), establishing the sought O(ε 6 ) expansion of J(ε) requires a O(ε 4 ) expansion of v ε on Γ ε Taking the previously known behavior [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] into account, an asymptotic expression for the trace on Γ ε of v ε for small ε will be sought in the form

v ε (ξ) = V 1 ( ξ) ε + V 2 ( ξ) ε 2 + V 3 ( ξ) ε 3 /2 + V 4 ( ξ) ε 4 /6 + o(ε 4 ) (ξ ∈ Γ ε , ξ ∈ S ) (24)
in terms of functions V 1 , . . . , V 4 defined on S . The determination of V 1 , . . . , V 4 , which constitutes the main step towards establishing an explicit expression for the expansion of J(ε), is based on expanding about ε → 0 an integral equation formulation for the scattered field v ε . This task is addressed in the next section.

Expansion of scattered field on the obstacle

Integral equation formulation of the forward problem. Let the Green's function G(x, ξ) associated with a unit point source located at x ∈ Ω be defined by the boundary-value problem

     (∆ ξ + k 2 )G(x, ξ) + δ(ξ -x) = 0 (ξ ∈ Ω) H(x, ξ) = 0 (ξ ∈ S N ), G(x, ξ) = 0 (ξ ∈ S D ), (25)
(with

H(x, ξ) = ∇ ξ G(x, ξ) • n(ξ)).
The trace on Γ ε of the scattered field v ε , defined by ( 5), then solves the boundary integral equation

1 2 v ε (x) + Γε H(x, ξ)v ε (ξ) dΓ ξ + Γε G(x, ξ)p(ξ) dΓ ξ = 0 (x ∈ Γ ε ), (26)

Expansion of the integral equation

To study the asymptotic behaviour of integral equation ( 26) as ε → 0, it is useful to introduce further scaled geometric quantities:

x = a + εx , r = εr , r = εr (x, ξ ∈ Γ ε ; x, ξ ∈ S ) (27)
in addition to definition [START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF] of ξ, and to split the Green's function [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF] according to:

G(x, ξ) = G(x, ξ) + G C (x, ξ), ( 28
)
where G is the (singular) acoustic fundamental solution for the free space, so that

G(x, ξ) = 1 4πr e ikr , ∇ ξ G(x, ξ) = ikr -1 4πr 3 e ikr r, (29)
with r = ξx and r = |r|, and the complementary part G C is smooth at ξ = x.

Lemma 1. Using the ansatz (24) (with functions V 1 , . . . , V 4 to be determined later) for the scattered field v ε , integral equation ( 26) has the following expansion about ε = 0:

4 a=1 ε a (a -1)! LV a (x) -F a (x) + o(ε 4 ) = 0 (x ∈ S ) ( 30
)
where L is the governing integral operator of exterior Neumann problems for the Laplace equation in the normalized domain R 3 \ B, defined by

f ∈ H 1/2 (S ) → Lf (x) = 1 2 f (x) - S 1 4πr 3 [r•n]f ( ξ) d Γξ (x ∈ S ), (31)
and F 1 (x), . . . , F 4 (x) are defined by

F 1 (x) = ∇u(a)• Y 1 (x) (32a) F 2 (x) = ∇ 2 u(a) : Y 2 (x) (32b) F 3 (x) = ∇ 3 u(a) :•Y 3 (x) -k 2 ∇u(a)• I : Y 3 (x) + x2 Y 1 (x) -2x• Y 2 (x) + 2F 3 (a) - ik 3 2π |B| u(a) + S k 2 4πr [r•n V 1 d Γξ (32c) F 4 (x) = ∇ 4 u(a) :: Y 4 (x) -3k 2 ∇ 2 u(a)• I : Y 4 (x) + x2 Y 2 (x) -2x• Y 3 (x) + 6F 4 (a) + 6∇F 3 (a) + ik 3 2π |B| ∇u(a) • x + S r•n 3k 2 4πr V 2 + ik 3 2π V 1 d Γξ (32d)
where ∇ k u(a) denotes the k-th order gradient of u evaluated at ξ = a, and having set

Y 1 (x) = - S 1 4πr n d Γξ Y 3 (x) = - S 1 4πr ξ ⊗ ξ ⊗ n d Γξ Y 2 (x) = - S 1 4πr ξ ⊗ n d Γξ Y 4 (x) = - S 1 4πr ξ ⊗ ξ ⊗ ξ ⊗ n d Γξ (33)
and

F 3 (z) = |B| ∇u(a) - S V 1 n d Γξ •∇ ξ G C (z, a) -|B| k 2 u(a)G C (z, a) (34a)
F 4 (z) = - S V 2 n d Γξ •∇ ξ G C (z, a) - S V 1 [ ξ ⊗ n] d Γξ : ∇ 2 ξ G C (z, a) (34b)
Proof. The proof rests on splitting the Green's function according to [START_REF] Kleiber | Parameter Sensitivity in Nonlinear Mechanics: Theory and Finite Element Computations[END_REF] in integral equation [START_REF] Jackowska-Strumillo | Topological optimization and inverse problems[END_REF], and invoking expansion

G(x, ξ) = 1 4πr ε -1 + ik 4π - k 2 r 8π ε - ik 3 r2 24π ε 2 + o(ε 2 ) (x, ξ ∈ Γ ε ), (35a)
H(x, ξ) = -[r•n 1 4πr 3 ε -2 + k 2 8πr + ik 3 12π ε + o(ε 2 ) (x, ξ ∈ Γ ε ). ( 35b
)
about ε = 0 of the singular free-space fundamental solution resulting from a substitution of the scaled distance function [START_REF] Kirsch | The domain derivative and two applications in inverse scattering theory[END_REF] into definition [START_REF] Kress | A hybrid method for two-dimensional crack reconstruction[END_REF]. First, one finds

1 2 v ε (x) + Γε H(x, ξ)v ε (ξ) dΓ ξ = 4 a=1 ε a (a -1)! LV a (x) - ε 3 2 S r•n k 2 4πr V 1 d Γξ - ε 4 6 S r•n 3k 2 4πr V 2 + ik 3 2π V 1 d Γξ + o(ε 4 ) (x ∈ S ) (36)
with the help of (21b), ansatz [START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF] and expansion (35b).

Second, combining the O(ε 4 ) Taylor expansion of p(ξ) = n( ξ) • ∇u(a + ε ξ) with expansion (35a) yields

G(x, ξ)p(ξ) = n( ξ)• ε -1 4πr ∇ + 1 4πr ξ•∇ 2 + ε 8π 1 r ( ξ ⊗ ξ) : ∇ 3 + 2ik ξ•∇ 2 -k 2 r∇ + ε 2 24π 1 r ( ξ ⊗ ξ ⊗ ξ) : ∇ 4 + 3ik( ξ ⊗ ξ) : ∇ 3 -3k 2 r ξ•∇ 2 -ik 3 r2 ∇ u(a) + o(ε 4 ) (37)
This expansion is next integrated over Γ ε and expressed with integrals over S through (21b). Invoking definitions [START_REF] Nintcheu Fata | A computational basis for elastodynamic cavity identification in a semi-infinite solid[END_REF], rearranging the O(ε 3 ) and O(ε 4 ) contributions thus arising by noticing that r2 = x2 + ξ2 -2x• ξ and r = r2 /r, and using explicit expressions

S ∇ 2 u(a) : ξ ⊗ n d Γξ = -|B| ∇ 2 u(a) : I = k 2 |B| u(a), S r2 ∇u(a)•n d Γξ = S (x 2 + ξ2 -2x• ξ)∇u(a)•n d Γξ = 2 |B| ∇u(a)• x,
obtained with the help of identities (23a-c), one obtains:

Γε G(x, ξ)p(ξ) dΓ ξ = -εY 1 (x)•∇ + ε 2 Y 2 (x) : ∇ 2 + ε 3 2 Y 3 (x) :•∇ 3 + k 2 I : Y 3 (x) + x2 Y 1 (x) -2x• Y 2 (x) •∇ + ik 3 2π |B| + ε 4 6 Y 4 (x) ::∇ 4 + 3k 2 I : Y 4 (x) + x2 Y 2 (x) -2x•Y 3 (x) :∇ 2 - ik 3 2π |B| ∇ u(a) (38)
Third, the complementary kernel G C (x, ξ) being smooth when x = ξ, the following Taylor expansions about ε = 0 hold for any x, ξ ∈ S :

G C (x, ξ) = G C (a, a) + ε (x•∇ x + ξ•∇ ξ )G C (a, a) + o(ε), (39a)
∇ ξ G C (x, ξ) = ∇ ξ G C (a, a) + ε (x•∇ x + ξ•∇ ξ )∇ ξ G C (a, a) + o(ε). (39b)
Invoking scaling (21a,b) and ansatz [START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF], one then obtains

Γε H C (x, ξ)v ε (ξ) dΓ ξ = S V 1 n d Γξ • (ε 3 + ε 4 x•∇ x )∇ ξ G C (a, a) + ε 4 S V 1 [ ξ ⊗ n] d Γξ : ∇ 2 ξ G C (a, a) + ε 4 S V 2 n d Γξ •∇ ξ G C (a, a) + o(ε 4 ),
by virtue of expansion (39b) and with

H C (x, ξ) = ∇ ξ G C (x, ξ)•n(ξ), and
Γε G C (x, ξ)p(ξ) dΓ ξ = Bε k 2 G C (x, ξ)u(ξ) -∇ ξ G C (x, ξ)•∇u(ξ) dV ξ = |B| k 2 u(a) (ε 3 + ε 4 x•∇ x )G C (a, a) -∇u(a)• (ε 3 + ε 4 x•∇ x )∇ ξ G C (a, a) + o(ε 4 )
with the help of the divergence theorem, identity [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] and expansions [START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Lochpositionierungskriterien[END_REF], (39a). Consequently, using definitions (34a,b) of F 3 , F 4 , one finds

Γε H C (x, ξ)v ε (ξ) + G C (x, ξ)p(ξ) dΓ ξ = -ε 3 F 3 (a) -ε 4 F 4 (a) + x•∇F 3 (a) + o(ε 4 ). ( 40
)
Lemma 1 now follows from setting to zero the sum of ( 36), ( 38) and [START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF], and reordering contributions according to powers of ε.

Functions V 1 , . . . , V 4 are to be found from equations LV a -F a = 0 arising from [START_REF] Litman | Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set[END_REF].

As F 3 , F 4 depend on V 1 , V 2 , this requires two stages: (i) solve LV a -F a = 0 (a = 1, 2) for V 1 , V 2 , (ii) inject the result into F 3 , F 4 and solve LV a -F a = 0 (a = 3, 4) for V 3 , V 4 . 4.2. Representation formulae for V 1 , V 2 Lemma 2. The functions V 1 and V 2 featured in expansion (24) of v ε are given by V 1 ( ξ; a) = ∇u(a)• U 1 ( ξ), V 2 ( ξ; a) = ∇ 2 u(a) : U 2 ( ξ) ( ξ ∈ S ), ( 41
)
where the vector function U 1 and the second-order tensor function U 2 do not depend on a and solve the integral equations

L U 1 (x) = Y 1 (x), L U 2 (x) = Y 2 (x) (x ∈ S ), (42)
Proof. From expansion (30) of integral equation ( 26), the governing equations for

V 1 , V 2 are [ LV 1 ](x) = ∇u(a)• Y 1 (x), [ LV 2 ](x) = ∇u 2 (a) : Y 2 (x) ( ξ ∈ S ),
with Y 1 , Y 2 and operator L defined in Lemma 1. As their right-hand sides depend linearly on ∇u(a), ∇ 2 u(a), representations (41) and equations (42) follow at once.

     ∆ U 1 = 0 ∇ U 1 •n = -n | ξ| 2 U 1 = O(1) ,      ∆ U 2 = 0 ∇ U 2 •n = -ξ ⊗ n | ξ| 2 U 2 = O(1)
,

( ξ ∈ R 3 \ B) ( ξ ∈ S ) (| ξ| → ∞). (43)
Lemma 3. Functions U 1 and U 2 defined by Lemma 2 verify the following integral identities:

S U 1 ⊗ n d Γξ = S n ⊗ U 1 d Γξ, ( 44a
) S U 2 ⊗ n d Γξ = S ξ ⊗ n ⊗ U 1 d Γξ. ( 44b
) S U 2 ⊗ ξ ⊗ n d Γξ = S ξ ⊗ n ⊗ U 2 d Γξ. ( 44c
)
Proof. The following relationship holds true for any pair (w , w ) of (scalar or tensor) functions that are harmonic in R 3 \ B and decay sufficiently fast at infinity:

S w (∇w •n) -w (∇w •n) d Γξ = 0. ( 45
)
Identities (44a,b,c) result from respectively applying (45) to pairs (w , w

) = ( U 1 , U 1 ), (w , w ) = ( U 2 , U 1 ) and (w , w ) = ( U 2 , U 2 )
, and invoking the boundary conditions (43) satisfied by U 1 , U 2 .

Representation formulae for

V 3 , V 4
Lemma 4. The functions V 3 and V 4 featured in expansion (24) of v ε are given by

V 3 ( ξ; a) = ∇ 3 u(a) :• U 3 ( ξ) + ∇u(a)•V 1 ( ξ) + α(a) (46a) V 4 ( ξ; a) = ∇ 4 u(a) :: U 4 ( ξ) + k 2 ∇ 2 u(a) : V 2 ( ξ) + γ(a) + β(a)• x (46b)
where the tensor functions U 3 , U 4 (respectively of order 3 and 4) do not depend on a and solve the integral equations

L U 3 (x) = Y 3 (x), L U 4 (x) = Y 4 (x) (x ∈ S ), ( 47
)
with Y 3 , Y 4 defined by [START_REF] Nintcheu Fata | A computational basis for elastodynamic cavity identification in a semi-infinite solid[END_REF], and the vector function V 1 and the second-order tensor function V 2 solve the integral equations

LV 1 (x) = Z 1 (x), LV 2 (x) = Z 2 (x) (x ∈ S ), (48)
having set

Z 1 (x) = 2x• Y 2 (x) -I : Y 3 (x) -x2 Y 1 (x) + S 1 4πr [ r•n ] U 1 d Γξ, Z 2 (x) = 3 2x• Y 3 (x) -I : Y 4 (x) -x2 Y 2 (x) + 3 S 1 4πr [ r•n ] U 2 d Γξ. ( 49
)
with Y 1 , Y 2 , Y 3 defined by [START_REF] Nintcheu Fata | A computational basis for elastodynamic cavity identification in a semi-infinite solid[END_REF]. Moreover, the scalars α(a), γ(a) and the vector β(a) featured in (46a,b) depend only on the sampling point a and are given by

α(a) = 2∇u(a)•A 11 •∇ ξ G C (a, a) -k 2 |B| u(a) 2G C (a, a) + ik 2π (50a) γ(a) = 6∇ 2 u(a) : A 21 •∇ ξ G C (a, a) + 6∇ 2 ξ G C (a, a) - ik 3 2π I : A 21 •∇u(a) (50b) β(a) = ∇u(a)•A 11 • 6∇ x ∇ ξ G C (a, a) + ik 3 2π I -6k 2 |B| u(a)∇ x G C (a, a) (50c)
with the constant tensors A 11 , A 21 given by (a)

A 11 = |B| I - S n ⊗ U 1 d Γξ, (b) A 21 = - S ξ ⊗ n ⊗ U 1 d Γξ. (51)
Proof. Upon substitution of expressions (41) of V 1 , V 2 into (32c) and (32d), and invoking definitions (49) of Z 1 , Z 2 and (51) of A 11 , A 21 , one obtains

F 3 (x) = ∇ 3 u(a) :•Y 3 (x) + k 2 ∇u(a)• Z 1 (x) + α(a), (52a)
F 4 (x) = ∇ 4 u(a) :: Y 4 (x) + k 2 ∇ 2 u(a) : Z 2 (x) + γ(a) + β(a)• x, (52b)
with α(a), γ(a), β(a) defined by (50a-c) and having exploited identity 54), together with definitions (47) of U 3 , U 4 and (48) of V 1 , V 2 , in (52a,b), the governing integral equations for V 3 , V 4 arising from expansion (30) of integral equation ( 26) take the form

|B| x + S r•n U 1 d Γξ = x•A 11 -I : A 21 , (53)
[ LV 3 ](x) = ∇ 3 u(a) :•[ L U 3 ](x) + k 2 ∇u(a) :•[ LV 1 ](x) + α(a) L(1) (x) (55a) [ LV 4 ](x) = ∇ 4 u(a) :: [ L U 4 ](x) + k 2 ∇ 2 u(a) : [ LV 2 ](x) + γ(a) L(1) (x) + β(a)• L ξ + U 1 (x) (55b)
Representations (46a), (46b) then follow directly from (55a) and (55b) by virtue of the fact that integral operator L is invertible.

Lemma 5. Functions U 3 and U 4 defined in Lemma 4 verify the integral identities:

S U 3 ( ξ) ⊗ n d Γξ = S ξ ⊗ ξ ⊗ n ⊗ U 1 d Γξ, ( 56a
) S U 4 ( ξ) ⊗ n d Γξ = S ξ ⊗ ξ ⊗ ξ ⊗ n ⊗ U 1 d Γξ, ( 56b
) S U 3 ( ξ) ⊗ ξ ⊗ n d Γξ = S ξ ⊗ ξ ⊗ n ⊗ U 2 d Γξ, (56c)
Proof. Properties (56a,b,c) result from the boundary conditions (57) satisfied by U 3 , U 4 and reciprocity identity (45), here applied to pairs (w , w ) = ( U 3 , U 1 ), (w , w ) = ( U 4 , U 1 ) and (w , w ) = ( U 3 , U 2 ), respectively.

Remark 2. Integral equations (47) correspond to the Laplace exterior Neumann problems on the normalized domain R 3 \ B defined by

     ∆ U 3 = 0 ∇ U 3 •n = -ξ ⊗ ξ ⊗ n | ξ| 2 U 3 = O(1) ,      ∆ U 4 = 0 ∇ U 4 •n = -ξ ⊗ ξ ⊗ ξ ⊗ n | ξ| 2 U 4 = O(1)
,

( ξ ∈ R 3 \ B) ( ξ ∈ S ) (| ξ| → ∞) (57)
Remark 3. Like U 1 and U 2 , the auxiliary tensor functions U 3 , U 4 , V 1 , V 2 are found by solving once and for all a set of integral equations, namely equations ( 47) and (48).

O(ε 6

) expansion of the misfit function

Arbitrarily shaped small scatterer

Building on the preliminary results established thus far, the main result of this article, namely the desired O(ε 6 ) expansion of J(ε), is now formulated.

Proposition 1. For a sound-hard scatterer represented by [START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF], i.e. of shape B and characteristic size ε, embedded in the acoustic reference medium Ω at a chosen location a in such a way that [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] holds, the O(ε 6 ) expansion of any objective function J(ε)

of format [START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF] with densities ϕ N (w R , w I , ξ) and ϕ D (w R , w I , ξ) twice differentiable w.r.t. their first two arguments is

J(ε; a) = J 6 (ε; a) + o(ε 6 ) ( 58
)
in terms of the sixth-order polynomial approximation

J 6 (ε; a) ≡ J(0) + T 3 (a)ε 3 + T 4 (a)ε 4 + T 5 (a)ε 5 + T 6 (a)ε 6 , (59)
with the coefficients T 3 (a), T 4 (a), T 5 (a) and T 6 (a) given by

T 3 (a) = Re ∇u•A 11 •∇û -|B| k 2 uû (a), (60a)
T 4 (a) = Re ∇ 2 û : A 21 •∇u + ∇ 2 u : A 21 •∇û (a), (60b)
T 5 (a) = Re 1 2 I 2 : ∇ 2 ∇u•∇û -k 2 uû + ∇û•B 11 •∇u + 1 2 ∇ 3 u :•A 31 •∇û + 2∇ 2 u : A 22 : ∇ 2 û + ∇ 3 û :•A 31 •∇u (a), (60c)
T 6 (a) = Re 1 6 I 3 : ∇ 3 ∇u•∇û -k 2 uû + 1 6 (∇ 4 u :: A 41 •∇û + ∇ 4 û :: A 41 •∇u) + 1 2 (∇ 3 u :•A 32 : ∇ 2 û + ∇ 3 û :•A 32 : ∇ 2 u) + k 2 ∇u•B 12 : ∇ 2 û + ∇ 2 u : B 21 •∇û - 1 2 |B| αû + 1 6 β•A 11 •∇û (a) + 1 2 a,b=R,I S N ϕ N,ab W a W b dΓ + S D ϕ D,ab Q a Q b dΓ . (60d)
In (60a-d), α(a), β(a) are defined by (50a,b), the function W is given by

W (x) = ∇ ξ G(x, a)•A 11 •∇u(a) -|B| k 2 G(x, a)u(a), (61)
and Q = ∇W •n, the constant tensors A 11 etc. are given by (51) and 42) and (48), and the constant tensors I 2 , I 3 are defined by

A 31 = - S ξ ⊗ ξ ⊗ n ⊗ U 1 d Γξ, A 41 = - S ξ ⊗ ξ ⊗ ξ ⊗ n ⊗ U 1 d Γξ, A 22 = - S ξ ⊗ n ⊗ U 2 d Γξ, A 32 = - S ξ ⊗ ξ ⊗ n ⊗ U 2 d Γξ, B 11 = - S V 1 ( ξ) ⊗ n d Γξ, B 21 = - S V 2 ( ξ) ⊗ n d Γξ, B 12 = - S V 1 ( ξ) ⊗ ξ ⊗ n d Γξ (62) in terms of the solutions U 1 , U 2 , V 1 , V 2 to equations (
I 2 = B ( ξ ⊗ ξ) d Vξ I 3 = B ( ξ ⊗ ξ ⊗ ξ) d Vξ. ( 63
)
Proof. The proof is straightforward, and consists in deriving an explicit form for expansion [START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF]. Each integral of ( 14) is expanded separately, as explained below. In particular, expansion of the second integral of ( 14) exploits the results of Section 4. Result (59), (60a-d) is finally obtained by adding all those contributions and rearranging the resulting expression.

(a) First integral of [START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF]. Using scaled coordinates according to [START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF] and (21a) and performing a straightforward Taylor expansion about ε = 0, one has

Bε ∇u•∇û -k 2 uû dV = ε 3 B ∇u•∇û -k 2 uû (a + ε ξ) d Vξ = ε 3 |B| + ε 5 2 I 2 : ∇ 2 + ε 6 6 I 3 :•∇ 3 ∇u•∇û -k 2 uû (a) + o(ε 6 ) (64)
with tensors I 2 (geometrical inertia of the normalized obstacle B) and I 3 given by ( 63).

(b) Second integral of [START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF]. Invoking ansatz [START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF], scaling (21b) and the counterpart for p of Taylor expansion [START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Lochpositionierungskriterien[END_REF], one has

Γε v ε p dΓ = ε 3 S V 1 ∇û(a)•n d Γξ + ε 4 S V 1 ξ•∇ 2 û(a) + V 2 ∇û(a) •n d Γξ + ε 5 2 S V 1 ξ ⊗ ξ : ∇ 3 û(a) + 2V 2 ξ•∇ 2 û(a) + V 3 ∇û(a) •n d Γξ + ε 6 6 S V 1 ξ ⊗ ξ ⊗ ξ : ∇ 4 û(a) + 3V 2 ξ ⊗ ξ : ∇ 3 û(a) + 3V 3 ξ•∇ 2 û(a) + V 4 ∇û(a) •n d Γξ + o(ε 6 ) (65)
Then, representations (41) of V 1 , V 2 and (46a,b) of V 3 , V 4 are inserted into (65), and integrations over S in the resulting formula are expressed in terms of the various constant tensors defined by ( 51) and (62), using when necessary relationships (44a,b,c) and (56a,b,c), and by invoking identities [START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF].

(c) Third and fourth integrals of [START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF]. The scattered field v ε at any point x ∈ S of the observation surface is given by the integral representation formula

v ε (x) = - Γε H(x, ξ)v ε (ξ) dΓ ξ - Γε G(x, ξ)p(ξ) dΓ ξ (x ∈ Ω ) (66)
Since G(x, ξ) and H(x, ξ) are smooth functions for any x ∈ Γ ε , the leading contribution to v ε (x) as ε → 0 results from a derivation formally identical to that of expansion [START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF], where (i) only the leading O(ε 3 ) order is retained, (ii) the complementary Green's function G C is replaced with the complete Green's function G, and (iii) the tensor A 11 is introduced. One hence obtains

v ε (x) = ε 3 W (x) + o(ε 3 ), q ε (x) = ε 3 Q(x) + o(ε 3 ) (x ∈ S)
i.e. [START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized Newton method[END_REF], with the function W given by (61) and Q = ∇W •n.

Remark 4. Integral identities (56a-d) render actual computation of U 3 , U 4 unnecessary, all the constant tensors needed in (59) being expressed in terms of U 1 , U 2 and V 1 , V 2 only.

Remark 5. The coefficient T 3 (a) associated with the leading O(ε 3 ) contribution to J(ε) is, as expected, the previously known topological derivative of J, i.e. ( 16).

Centrally-symmetric scatterer

When B has central symmetry (i.e. is such that ξ ∈ B ⇔ -ξ ∈ B), the constant tensors I 3 defined by (63), A 21 defined by (51b) and A 41 , A 32 , B 21 , B 12 defined by (62) vanish, as shown in Appendix A. Moreover, actual computation of the auxiliary solution V 2 is no longer needed for setting up expansion (59). Consequently, Proposition 1 becomes:

Proposition 2. When the obstacle (20) has a centrally-symmetric shape B (i.e. ξ ∈ B ⇔ -ξ ∈ B), expansion (59) holds with coefficients T 3 , T 5 still given by (60a,c) and

T 4 (a) = 0, (67a)
T 6 (a) = Re - 1 2 |B| α(a)û(a) + 1 6 β(a)•A 11 •∇û(a) + 1 2 a,b=R,I S N ϕ N,ab W a W b dΓ + S D ϕ D,ab Q a Q b dΓ . (67b)
Many simple obstacle shapes, e.g. spheres, ellipsoids or rectangular boxes, are centrally-symmetric, which makes this simplification quite useful.

Spherical scatterer

The special case of a spherical sound-hard obstacle B ε (where B is the unit ball, S the unit sphere and |B| = 4π/3), which exhibits central symmetry and also permits further analytical treatment, is now considered.

The constant tensor I 2 defined by (63) is easily found to be given by

I 2 = 4π 9 I (68)
Moreover, the integral equations for the auxiliary unknowns U 1 , U 2 , V 1 are solvable in closed form. First, U 1 , U 2 are found, by applying the method of separation of variables in spherical coordinates to the exterior Laplace problems (43), to be given by

U 1 (x) = 1 2 x, U 2 (x) = 1 3 (x ⊗ x) + 2 9 I (69)
Then, inserting U 1 given by (69) into the right-hand side of integral equation (48a) and noting that n = -ξ on the sphere S , one finds after some manipulation that

V 1 satisfies equation LV 1 (x) = - 5 16π S x -ξ ξ d Γξ = 1 3 x ( 70
)
and is hence given, with the help of identity (54), by

V 1 ( ξ) = 1 3 ( U 1 ( ξ) + ξ) = 1 2 ξ ( 71
)
Based on solutions (69) and (71), the various constant tensors featured in formulae (60a,c) and (67b) are now obtained, using elementary integration methods, as

A 11 = 2πI , B 11 = π 3 I , A 22 = 4π 45 J + 8π 27 (I ⊗ I) , A 31 = 2π 15 J , (72)
where the fourth-order tensor J is defined by

J ijk = δ ij δ k + δ ik δ j + δ i δ jk (73)
On substituting these values into (60a,c) and (67b) and recalling result (67a), the O(ε 6 ) expansion of J(ε) is hence given a more explicit form:

Proposition 3. When the obstacle of Proposition 1 is spherical (with B chosen as the unit sphere), expansion (59) holds, with coefficients T 3 , . . . , T 6 now given by

T 3 (a) = Re 2π∇u•∇û - 4π 3 k 2 uû (a), (74a)
T 4 (a) = 0, (74b)

T 5 (a) = Re 4π 9 ∇ 2 : ∇ 2 û - 3π 5 k 2 ∇u•∇û + 88π 135 k 4 uû (a), (74c)
T 6 (a) = Re - 2π 3 α(a)û(a) + π 3 β(a)•∇û(a) + 1 2 a,b=R,I S N ϕ N,ab W a W b dΓ + S D ϕ D,ab Q a Q b dΓ . (74d)

Extension to several scatterers

Expressions (60a-d) of T 3 (a), . . . , T 6 (a) are predicated on the assumption of a single impenetrable scatterer characterized by its shape B, size ε and location a. However, this result can be extended to the case of K > 1 scatterers B (k) ε defined according to

B (k) ε (a (k) ) = a (k) + εB (k) (1 ≤ k ≤ K) ( 75
)
where a (k) and B (k) are the centre and (normalized) shape of the k-th scatterer, and the size parameter ε is the same for all K scatterers. To help present this generalization in a compact way, the following notational convention will be used: a superscript '(k)' appended to any previously defined symbol (e.g.

V (k) 1 , A (k) 11 , I (k)
3 ) will refer to quantities associated with the previous single-scatterer analysis, with B ε replaced by B (k) ε . Proposition 4. For a set of K sound-hard obstacles of form (75) embedded in the reference medium Ω at prescribed locations a (1) , . . . , a (K) , let J(ε; a (1) , . . . , a (K) ) be defined by [START_REF] Céa | The shape and topological optimization connection[END_REF], with 1) , . . . , a (K) ) denoting the scattered field induced by the K objects. Densities ϕ N (w R , w I , ξ), ϕ D (w R , w I , ξ) are assumed to be twice differentiable w.r.t. their first two arguments. The O(ε 6 ) expansion of J(ε) is J(ε; a (1) , . . . , a (K) ) = J(0)

Ω ε ≡ Ω\ B(1) ε ∪. . .∪ B(K) ε and v ε ≡ v ε (ξ; a (
+ K k=1 ε 3 T (k) 3 (a (k) ) + ε 4 T (k) 4 (a (k) ) + ε 5 T (k) 5 (a (k) ) + ε 6 T (k) 6 (a (1) , . . . , a (K) ) + o(ε 6 ) (76)
where

T (k) 3 , T (k) 4 , T (k) 5
are given by (60a-c) with shape B = B (k) , and

T (k) 6
is given by

T (k) 6 (a (1) , . . . , a (K) ) = T (k) 6 (a (k) ) + =k Re ∇E ( ) 3 •A k 11 •∇û -k 2 |B (k) |E ( ) 3 û (a (k) ) + 1 2 =k a,b=R,I S N ϕ N,ab W ( ) a W (k) b dΓ + S D ϕ D,ab Q ( ) a Q (k) b dΓ ( 77
)
where W ( ) are defined by ( 61) with a = a ( ) and B = B ( ) , and with E ( )

3 , E ( ) [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF] given by (83a,b).

Proof. The O(ε 6 ) expansion of J(ε) is sought on the basis of

J(ε) -J(0) = K k=1 B (k) ε ∇u•∇û -k 2 uû dV - K k=1 Γ (k) ε v ε p dΓ + 1 2 a,b=R,I S N ϕ N,ab v ε a v ε b dΓ + S D ϕ D,ab q ε a q ε b dΓ + o(|v ε | 2 L 2 (S N ) , |q ε | 2 L 2 (S D ) ). ( 78
)
The first term in the right-hand side of ( 78) is clearly a sum of contributions of the form (64) arising from each scatterer. Moreover, on noting that the integral representation (66) is a sum of integrals over each scatterer and revisiting the analysis of Section 4, the leading O(ε 3 ) contribution to v ε is simply the corresponding sum of contributions [START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized Newton method[END_REF], i.e.:

v ε (ξ) = ε 3 K k=1 W (k) (ξ) + o(ε 3 ) , q ε (ξ) = ε 3 K k=1 Q (k) (ξ) + o(ε 3 ) (ξ ∈ S) (79)
where W (k) is defined by (61) with a = a (k) and B = B (k) , and with

Q (k) = ∇W (k) •n.
The leading contribution of the last two integrals of (78), of order O(ε 6 ), then stems directly from (79). Finally, to evaluate the second integral of (78), an expansion of v ε on each scatterer, of the form

v ε (ξ) = ε V (k) 1 ( ξ)+ε 2 V (k) 2 ( ξ)+ε 3 V (k) 3 ( ξ)+ε 4 V (k) 4 ( ξ)+o(ε 4 ) (ξ ∈ Γ (k) ε , ξ ∈ S (k) ) (80) is again postulated. It is expected that ( V (k) 1 , . . . , V (k) 4 ) = (V (k) 1 , . . . , V (k)
4 ) due to multiple-reflection effects between scatterers. The governing integral equation for v ε is ( 26) with all integrals over Γ ε changed to sums of integrals over the Γ

(k) ε , i.e. 1 2 v ε (x) + Γ (k) ε H(x, ξ)v ε (ξ) + G(x, ξ)p(ξ) dΓ ξ + =k Γ ( ) ε H(x, ξ)v ε (ξ) + G(x, ξ)p(ξ) dΓ ξ = 0 (x ∈ Γ (k) ε , 1 ≤ k ≤ K). (81) The V (k) 1 , . . . , V (k) 4
are to be found by inserting (80) into (81) and expanding the resulting equations in powers of ε. A comparison with [START_REF] Jackowska-Strumillo | Topological optimization and inverse problems[END_REF] shows that the first line in (81) in terms of the Green's function G defined by [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF]. As a result, the O(ε 6 ) expansion of J(ε) established in Section 5 is almost completely explicit, the only non-explicit components being the auxiliary solutions U 1 , V 1 etc., which must be computed numerically except for simple shapes of the hypothetical scatterer B such as the spherical shape discussed in Section 5.3.

In practice, this explicit character is retained only for geometrically simple configurations (Ω, S N , S D ) such that the corresponding Green's function is known analytically. For instance, for the acoustic half-space Ω = {ξ | ξ 3 ≤ 0} bounded by

S = {ξ | ξ 3 = 0}, it is well-known that G C (x, ξ) = ± 1 4πr e ikr , with r = |ξ -x| , x = (x 1 , x 2 , -x 3 ) (89)
where the '+' and '-' sign correspond to the cases S N = S, S D = ∅ (Neumann) and

S D = S, S N = ∅ (Dirichlet).
For configurations where the Green's function is not known, the free and adjoint fields may be computed by solving the boundary integral equations [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF][START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF]

L(u, p) (x) = F(u D , p D ) (x), (90)
L(û, p) (x) = F(-ϕ D,q , ϕ N,u ) (x), (91)
with x ∈ S and where the linear integral operator L(f, g) and the right-hand side functional F(f D , g D ) are defined by

L(f, g) (x) = 1 2 f (x) + S N H(x, ξ)f (ξ) dΓ ξ - S D G(x, ξ)g(ξ) dΓ ξ (x ∈ S), (92a)
F(f D , g D ) (x) = S D H(x, ξ)f D (ξ) dΓ ξ - S N G(x, ξ)g D (ξ) dΓ ξ (x ∈ S), (92b)
and invoking subsequent integral representation formulae. Moreover, the pair (W, Q) associated with the leading O(ε 3 ) contribution of (v ε , q ε ) on S, defined by (61), and the complementary kernel pair G C (z, ξ), H C (z, ξ) , defined by [START_REF] Kleiber | Parameter Sensitivity in Nonlinear Mechanics: Theory and Finite Element Computations[END_REF] and featured in T 6 through α, β given by (50a,b), are respectively governed by integral equations

L(W, Q) (x) = k 2 |B| u(a)G(x, a) -∇u(a)•A 11 •∇G(x, a) (x ∈ S) (93)
L G C (z, •), H C (z, •) (x) = -F G(z, •), H(z, •) (x) (x ∈ S, z ∈ Ω) (94)
Topological sensitivity and low-frequency asymptotics. In cases featuring only one characteristic length (one scatterer embedded in an unbounded medium and illuminated by an incident plane wave), the present approach based on topological sensitivity is essentially similar to low-frequency direct and inverse scattering problems [START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF][START_REF] Dassios | Low frequency scattering[END_REF] with a cost function based on measurements taken at infinity. Such situations correspond in this paper to considering one single scatterer and setting the complementary Green's function G C to zero. Both approaches can then be reformulated in terms of expansions with respect to kε 1 (with k fixed and ε small here, but k small and ε fixed in [START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF][START_REF] Dassios | Low frequency scattering[END_REF]).

For example, results of Section 4 for the asymptotic behavior of v ε on Γ ε can, in the case G C = 0, be recovered (after rescaling and, for V 3 and V 3 , some manipulation) from the low-frequency analysis of the Neumann acoustic problem given in Sec. 3.A of [START_REF] Dassios | Low frequency scattering[END_REF]. However, if other characteristic lengths (e.g. the distance of a scatterer to a free surface, the size of a bounded region, or the radius of curvature of wave fronts) are present in addition to the vanishing size of the scatterer, the equivalence (up to scaling) between the topological-sensitivity and low-frequency approaches is no longer true. Direct vs. adjoint approaches for topological sensitivity. Topological sensitivity has formal similarities with the more traditional areas of parameter sensitivity [START_REF] Kleiber | Parameter Sensitivity in Nonlinear Mechanics: Theory and Finite Element Computations[END_REF] or shape sensitivity [START_REF] Sokolowski | Introduction to shape optimization[END_REF]. Like first-order parameter or shape sensitivity fomulae, the topological derivative T 3 associated with the leading O(ε 3 ) contribution to J(ε) is expressed as a bilinear combination of the free and adjoint fields, see also e.g. [START_REF] Bonnet | Topological sensitivity for 3D elastodynamic and acoustic inverse scattering in the time domain[END_REF][START_REF] Céa | The shape and topological optimization connection[END_REF][START_REF] Dominguez | Time domain topological gradient and time reversal analogy: an inverse method for ultrasonic target detection[END_REF][START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF][START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF]. Moreover, setting up the O(ε 6 ) expansion of J(ε), and particularly the highest-order coefficient T 6 , requires the 'direct topological sensitivities' W, Q associated with the leading O(ε 3 ) contributions to v ε and q ε on the measurement areas, in addition to the free and adjoint fields. This is reminiscent of the fact that second-order parameter or shape sensitivity fomulae can be cast as bilinear combinations of the free and adjoint fields and their first-order sensitivities. It is nevertheless important to keep in mind that topological and shape sensitivities are related but distinct concepts, as emphasized in [START_REF] Céa | The shape and topological optimization connection[END_REF].

Here, it would have been possible to establish the O(ε 6 ) expansion of J(ε) on the basis of (9) rather than [START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF], i.e. without recourse to the adjoint solution [START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF]. This alternative 'direct' approach requires O(ε 6 ) expansions of v ε on S N and q ε on S D , i.e. the actual computation of auxiliary solutions W 4 , W 5 , W 6 in addition to W 3 = W defined in [START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized Newton method[END_REF], which can be obtained by expanding integral representation (66) to order O(ε 6 ). Such high-order expansions of the field quantities are given, to arbitrary order and for various physical contexts, by Ammari and Kang [START_REF] Ammari | Reconstruction of small inhomogeneities from boundary measurements[END_REF].

Summation over experiments. The results of Sections 3 and 6 are easily generalized to cost functions J defined on the basis of N experiments, by setting

J (Ω ) = N e=1 J e (Ω )
where for each e J e (Ω ) is a cost function of format (6) associated with the e-th experiment. Each experiment gives rise to a free field u e and an adjoint field ûe . Hence, expansion (59) of J(ε) follows by summing over e the expansions (59) for each J e .

Numerical implementation

A simple approximate global search procedure

Expansions of the form (59) offer the option of minimizing the polynomial approximation J 6 (ε; a) of J(ε; a) for sampling points a chosen a priori. This task is, for each sampling point, simple and computationally very light. It can therefore be performed for locations a spanning a fine search grid G, thereby defining an approximate global search procedure over the spatial region sampled using G. The best estimate of the unknown scatterer B true yielded by this procedure is defined by the location a = x est and size ε = R est achieving the lowest value of J 6 (ε; a) over G, i.e. given by

x est = arg min a∈G Ĵ6 (a), R est = R(x est ), (95)

Numerical results

To demonstrate the proposed approximate search procedure and thereby demonstrate the usefulness of the O(ε 6 ) expansion of J(ε; a), the identification of an impenetrable object embedded in an acoustic medium occupying the half-space Ω = {ξ | ξ 3 ≤ 0} is considered, as depicted on Fig. 1. A homogeneous Neumann condition is assumed on the surface S = S N = {ξ | ξ 3 = 0}. Under these conditions, the relevant Green's function G is, as mentioned earlier, explicitly known and given by ( 28), ( 29) and (89).

where the overbar indicates complex conjugation. The synthetic data u obs e are computed by means of a direct boundary element method (BEM) wherein the boundary of B true is meshed using 600 eight-noded boundary elements.

The scatterer B true to be identified is centered at x true = (2.05a, 1.25a, -3.05a). Three geometries are considered for B true (Figure 1): a sphere (S) of radius 0.5a, a horizontally elongated ellipsoid (E) with semiaxes (a, 0.5a, 0.5a), and a banana-shaped scatterer (B) obtained by applying the transformation ξ 3 -→ ξ 3 -(ξ 1 -x true 1 ) 2 /2a to the ellipsoid of semiaxes (a, 0.25a, 0.25a), the semiaxes being aligned in both cases with the (ξ 1 , ξ 2 , ξ 3 ) coordinates. For comparison purposes, the 'true' radius R true is defined as the radius of the sphere having the same volume as B true , i.e. R true = 2 -1 a, 2 -2/3 a, 2 -4/3 a, respectively, for (S), (E) and (B). Three wavenumbers ka = 0.5, 1, 2 have been considered. Note that x true ∈ G: the sampling point a ∈ G closest to x true has coordinates (2a, 1.2a, -3.2a) and is separated from x true by a distance √ 11a/20 ≈ .166a.

Determination of obstacle size (known location, noise-free data). In a preliminary numerical experiment, the size of an obstacle of known location is estimated via the computation of R(x true ) defined by (96a), where J 6 (ε; x true ) is defined in terms of a trial spherical scatterer, i.e. using coefficients T 3 (x true ), . . . , T 6 (x true ) given by (74a-d).

Results for the relative error R(x true )/R true -1 on the obstacle size estimation obtained using noise-free synthetic data are given in Table 1 for all of the previously-defined obstacle configurations, testing configurations and wavenumbers. These results indicate in particular that the size estimation accuracy decreases as the frequency increases, and is relatively insensitive to the density of the testing and measurement grids.

Approximate global search procedure, noise-free data. The approximate global search procedure defined in Section 8.1 has been performed on a search grid G of 51×51×25 = 65025 regularly spaced sampling points spanning the 3-D box-shaped region defined by -10a ≤ x 1 , x 2 ≤ 10a, -10a ≤ x 3 ≤ -0.4a. The polynomial approximation J 6 (ε; a) of cost function (97) associated with a trial spherical scatterer, i.e. with coefficients T 3 (a), . . . , T 6 (a) again given by (74a-d), has been set up for all 65025 sampling points a ∈ G of the search grid thus defined and using the explicit Green's function.

The obstacle radius estimation R est defined by (95) obtained for all of the previously-defined obstacle configurations, testing configurations and wavenumbers are compared to R true in Table 2, using noise-free synthetic data. The lower-frequency case k = 0.5a is again seen to yield the most accurate estimation of R true . For all cases displayed in Table 2, the identified obstacle location a est is the grid point closest to x true , i.e. a estx true = √ 11a/20 ≈ .166a. A comparison of the results of Tables 1 and2 shows that, not surprisingly, R(x true ) is usually a slightly more accurate estimation of R true than R est . For cases (E) and (B), featuring 'true' scatterer shapes that increasingly deviate from the trial spherical shape, the accuracy for the 'equivalent radius' R true is nonetheless similar to that obtained for case (S). The size estimation accuracy is, again, seen to decrease as the frequency increases, and is relatively insensitive to the density of the testing and measurement grids. Even the 2 × 2 testing configuration, featuring only 4 sources and 9 sensors, yields good results when applied to error-free data.

For comparison purposes, the obstacle radius estimation R est has also been computed, using the same sampling grid G, for true scatterers (S'), (E') and (B') with the same location and shape as (S), (E) and (B) and size reduced by a factor 2.5 (e.g. (S') is a sphere of radius 0.2a), using testing configurations 5 × 5, 10 × 10 and wavenumbers ka = 0.5, 1, 2. The O(ε 6 ) expansion of J(ε), and hence the radius estimates, are expected to be more accurate for this set of smaller true scatterers. On comparing the relative error on R est obtained from these computations, presented in 2, results are seen to conform to this expectation, except in some of the cases with ka = 0.5, and otherwise follow the same mentioned trends.

In the results presented so far, only the minimum J min 6 = Ĵ6 (x est ) achieved by Ĵ6 (a) over a ∈ G was considered. However, another interesting outcome of the numerical experiments performed is that values of Ĵ6 (a) close to the minimum Ĵmin 6 are found to occur only at grid points close to a est and to yield optimal radii R(a) similar to R est . To illustrate this finding, iso-surfaces of Ĵ6 (a) for Ĵ6 = ζJ min 6 , with ζ = 0.6 , 0.7 , 0.8 , 0.9, depicted on Figure 2 for obstacle configuration (E) and testing configuration 20 × 20, are seen to shrink to a small neighbourhood of x true as Ĵ6 approaches J min 6 .

Approximate global search procedure, noisy data. Finally, the effect of data errors on the approximate search procedure is examined. The (synthetically) measured total field u obs where η e,m and η e,m are uniform random numbers with zero mean and 0.05 standard deviation. The measurement residuals u e (x m ) -u obs e (x m ) being on average of much smaller magnitude than the measured total field, especially at the lower frequency ka = 0.5, they are severely affected by the above-defined, relatively small, perturbation of the total field. Estimations R est and x est for true scatterer configurations (S), (E) and (B) have been computed, using the same sampling grid G, for testing configurations BH, 5 × 5, 10 × 10, 20 × 20 and wavenumbers ka = 0.5, 1, 2. The relative error on R est and the distance a estx true resulting from these computations are displayed in Table 4.

Table 4. Distance a estx true and relative error R est /R true -1 for obstacles (S), (E) and (B) of unknown location, testing configurations BH, 5×5, 10×10 and 20×20, and synthetic data with 5% noise on total field. Where a estx true is unacceptably large, relative error on obstacle size is deemed irrelevant and not shown. A comparison of these results with those for error-free data (Table 2) shows that mildly perturbed values of R est and x est are obtained for ka = 1, 2, whereas the deterioration of accuracy is much stronger for ka = 0.5 when the coarser testing configurations 5 × 5 and BH are used.

Conclusion

In this article, extending previous work on topological sensitivity, a methodology for expanding to order O(ε 6 ) a generic misfit cost function associated with the identification of obstacles of characteristic size ε has been developed. Although presented for the specific case of sound-hard obstacles in linear acoustic media, the approach is generic and is expected to yield similar expansions for other cases, e.g. penetrable obstacles in acoustic, elastic or electromagnetic. A non-iterative fast approximate global search strategy based on a very simple exploitation of the O(ε 6 ) expansion has been proposed and demonstrated on numerical experiments to correctly identify a single scatterer, even in the case of noisy data. Future work include the extension of this approach to penetrable obstacles and cracks, and to the case of multiple scatterers for which the present approximate global search algorithm cannot be applied without modification.

S

  ξ⊗n d Γξ = -|B| I,(c)S ξ⊗ ξ⊗n d Γξ = 0,[START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF]

Remark 1 .

 1 Integral equations (41) correspond to the following Laplace exterior Neumann problems on the normalized domain R 3 \ B (with symbols ∇, ∆ indicating differentiation with respect to normalized coordinates ξ):

  which holds by virtue of r = ξx and definitions (51), for equation (52b). Moreover, identities L(1) (x) = 1 , L ξ + U 1 (x) = x(54) hold true by virtue of the fact that functions f 0 ( ξ) = 1 and f 1 ( ξ) = ξ solve the integral equation associated with interior Neumann problems on domain B for the Laplace equation with respective boundary data p D 0 = 0 and p D 1 = -n. Invoking (

Figure 1 .

 1 Figure 1. Identification of a spherical or ellipsoidal sound-hard scatterer in a acoustic half-space: geometry and notation.

  with functions Ĵ6 (a) and R(a) defined through a partial minimization of J 6 (ε; a) w.r.t. ε, i.e.: Ĵ6 (a) = min ε J 6 (ε; a), R(a) = arg min ε J 6 (ε; a). (96)

Four synthetic testing configurations (labelled 2 × 2 , 5 × 5 , 10 ×e

 225510 10 and 20 × 20 in the sequel) are defined, where the square region {ξ | -5a ≤ ξ 1 , ξ 2 ≤ 5a, ξ 3 = 0} of S is divided into 2 × 2, 5 × 5, 10 × 10 and 20 × 20 squares, respectively (a being a reference length). Acoustic point sources x e and sensors x m are located at all centers and vertices, respectively, of the above-defined square grids, so that configurations 2×2, 5×5, 10×10 and 20×20 feature N = 4, 25, 100, 400 sources and M = 9, 36, 121, 441 sensors, respectively. A fifth testing configuration (reminiscent of borehole measurements and hereinafter labelled BH) consists of two sets of 10 vertically-aligned and evenly spaced point sources x e = 5a(-1, 1, -e/10) (1 ≤ e ≤ 10), x e = 5a(1, -1, 1 -e/10) (11 ≤ e ≤ 20) and two sets of 10 vertically-aligned and evenly spaced receiversx m = 5a(1, 1, -m/10) (1 ≤ m ≤ 10), x m = 5a(-1, -1, 1 -m/10) (11 ≤ m ≤[START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF], so that M = N = 20. All of these configurations define situations of limited aperture since data is not available in all directions around the sought scatterer.The N point sources are applied in sequence, thus defining a set of N synthetic experiments. The identification is formulated in terms of the least-squares cost function J (and u e denote the acoustic fields induced by point source x e for the 'true' and 'trial' configurations Ω true = Ω \ (B true ∪ Γ true ) and Ω = Ω \ (B ∪ Γ ). The free and adjoint fields associated with the e-th experiment are given by u e (ξ) = G(x e , ξ) ûe (ξ) = M m=1 u e (x m ) -u obs e (x m ) G(x m , ξ)

Figure 2 . 6 ,

 26 Figure 2. Iso-surfaces of Ĵ6 (a) for Ĵ6 = ζJ min 6 , with ζ = 0.6 (top left), ζ = 0.7 (top right), ζ = 0.8 (bottom left) and ζ = 0.9 (bottom right), for obstacle configuration (E), testing configuration 5 × 5 and noise-free data. The location of the true scatterer and of the iso-surfaces are emphasized via projections onto three orthogonal planes.

e

  in cost function (97) has been replaced with a perturbed version ũobs e such that Re ũobs e (x m ) = (1 + η e,m )Re u obs e (x m ) , Im ũobs e (x m ) = (1 + η e,m )Im u obs e (x m )

RR

  -01 -1.37e-01 -3.07e-01 10 × 10 -4.15e-02 -1.03e-01 -3.08e-01 20 × 20 -5.28e-02 -9.97e-02 -3.08e-01 (E) a estx true BH 8est /R true -1 BH ---8.35e-02 -3.31e-01 5 × 5 1.20e-02 -1.39e-01 -3.92e-01 10 × 10 -3.39e-02 -1.37e-01 -3.92e-01 20 × 20 -1.49e-03 -1.36e-01 -3.93e-01 est /R true -1 BH -----1.83e-01 5 × 5 ---2.44e-02 -2.22e-01 10 × 10 ---1.03e-02 -2.13e-01 20 × 20 ---3.69e-02 -2.11e-01

Table 1 .

 1 Relative error R(x true )/R true -1 for radius estimates R(x true ), for obstacles (S), (E) and (B) of known location, testing configurations BH, 5×5, 10×10 and 20×20, and noise-free synthetic data.

		ka = 0.5	ka = 1	ka = 2
	(S) BH	-2.73e-02 -9.70e-02 -2.75e-01
	5 × 5	-2.09e-02 -1.03e-01 -3.02e-01
	10 × 10 -2.07e-02 -1.03e-01 -3.03e-01
	20 × 20 -2.05e-02 -1.03e-01 -3.03e-01
	(E) BH	-3.26e-02 -1.33e-01 -3.53e-01
	5 × 5	-1.62e-02 -1.43e-01 -3.93e-01
	10 × 10 -1.56e-02 -1.42e-01 -3.93e-01
	20 × 20 -1.52e-02 -1.42e-01 -3.93e-01
	(B) BH	-4.31e-03 -4.76e-02 -2.03e-01
	5 × 5	1.39e-02 -4.80e-02 -2.11e-01
	10 × 10	1.46e-02 -4.73e-01 -2.10e-01
	20 × 20	1.50e-02 -4.69e-01 -2.10e-01

Table 2 .

 2 Relative error R est /R true -1 on radius estimation R est for obstacles (S), (E) and (B) of unknown location, testing configurations BH, 5×5, 10×10 and 20×20, and noise-free synthetic data. A distance a estx true = a √ 11/20 is found for all cases.

		ka = 0.5	ka = 1	ka = 2
	(S) BH	-2.19e-02 -9.10e-02 -2.77e-01
	2 × 2	-1.28e-02 -9.99e-02 -3.06e-01
	5 × 5	-1.02e-02 -9.88e-02 -3.08e-01
	10 × 10 -9.70e-03 -9.86e-02 -3.08e-01
	20 × 20 -9.45e-03 -9.85e-02 -3.08e-01
	(E) BH	-4.50e-02 -1.27e-01 -3.47e-01
	2 × 2	-8.31e-03 -1.39e-01 -3.92e-01
	5 × 5	-4.51e-03 -1.38e-01 -3.93e-01
	10 × 10 -3.58e-03 -1.37e-01 -3.93e-01
	20 × 20 -3.10e-03 -1.37e-01 -3.93e-01
	(B) BH	4.87e-03 -4.09e-02 -2.09e-01
	2 × 2	2.35e-02 -4.30e-02	-2.14e-01
	5 × 5	2.66e-02 -4.09e-02	-2.13e-01
	10 × 10	2.76e-02 -4.00e-02	-2.12e-01
	20 × 20	2.81e-02 -3.95e-02	-2.12e-01

Table 3 .

 3 Relative error R est /R

		ka = 0.5	ka = 1	ka = 2
	(S') 5 × 5	7.54e-03 -1.19e-02 -7.83e-02
	10 × 10	7.85e-03 -1.16e-02 -7.82e-02
	(E') 5 × 5	2.32e-02 -5.82e-03 -9.60e-02
	10 × 10	2.38e-02 -5.27e-03 -9.56e-02
	(B') 5 × 5	4.07e-02 -2.30e-02 -3.24e-02
	10 × 10	4.16e-02 -2.37e-02 -3.17e-02

true -1 on radius estimation R est for smaller obstacles (S'), (E') and (B') of unknown location, testing configurations 5 × 5 and 10 × 10, and noise-free synthetic data. A distance a estx true = a √ 11/20 is found for all cases.

Table 3 ,

 3 with corresponding results of Table

constitute the contribution arising due to scatterer B (k) ε in isolation. The expansion in ε of that contribution therefore coincides with that established in Section 3.3 for the single-scatterer case. Besides, the sum of integrals in the second line of (81), which synthesizes the influence of scatterers B ( )

ε , can readily be shown by means of a calculation similar to that leading to [START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF] to have the expansion

where the scalar functions

4 are defined for any x = a ( ) by E ( )

Since contributions (82) are of order O(ε 3 ), the O(ε) and O(ε 2 ) contributions to equation (81) are not affected by the scatterers B ( )

and one therefore has

Moreover, the form assumed by the supplementary contributions ( 82) is such that results of Section 4 still apply provided every occurrence of α(a), γ(a),

(a (k) )

(85) The supplementary terms (contributions of B ( ) ε , = k) are the only manifestations of multiple reflections between scatterers arising in this analysis. The auxiliary unknowns

are then given by (46a,b) with replacements (85), i.e. by

Proposition 4 then follows from collecting results (78), ( 79), ( 84), ( 85) and (86a,b) and revisiting the analysis of Sections 4 and 5.

Discussion

Computational issues. The free and adjoint fields are given by the following explicit formulae The mapping σ is in particular such that

Then, I 3 vanishes by virtue of

The other above-mentioned tensors also vanish, as a consequence of symmetry properties of the solutions to integral equations ( 42) and (48,b). Specifically, U 2 , V 2 are symmetric while U 1 , V 1 are skew-symmetric. For example, to prove that U 1 is skew-symmetric, let U even 1 and U odd 1 , the even and odd parts of U 1 , be defined by:

These definitions imply that

Now, on inserting the decomposition U 1 = U even 1 + U odd 1 in integral equation (47a), writing the resulting equations for a pair of symmetrical collocation points x and σ x (x ∈ S ), using property (A.3), and noting that

the following pair of integral equations is arrived at:

On taking the sum and the difference of equations (A.4), an equivalent system of integral equations (A.4) is obtained: Leven