Marc Bonnet
email: mbonnet@ensta.fr

INVERSE ACOUSTIC SCATTERING USING HIGH-ORDER SMALL-INCLUSION EXPANSION OF MISFIT FUNCTION

Keywords: 2010 Mathematics Subject Classification. Primary: 31A10, 31A25, 35J05, 58J37 Inverse scattering, Helmholtz equation, topological derivative, asymptotic expansion, volume integral equation

Inverse acoustic scattering using high-order small-inclusion expansion of misfit function

1. Introduction. Inverse scattering has been intensively investigated over the last quarter century, in particular due to the development of qualitative, sampling-based, methods [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF][START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF][START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF] that offer robust and computationally effective alternatives to more-classical techniques based on e.g. PDE-constrained minimization or successive linearizations. Linear sampling and factorization methods have in particular been extensively studied. In adddition, the concept of topological derivative has been investigated in a variety of inverse scattering situations towards defining imaging hidden objects, see e.g. [START_REF] Ammari | Stability and resolution analysis for a topological derivative based imaging functional[END_REF][START_REF] Bellis | A FEM-based topological sensitivity approach for fast qualitative identification of buried cavities from elastodynamic overdetermined boundary data[END_REF][START_REF] Bellis | Acoustic inverse scattering using topological derivative of far-field measurements-based L 2 cost functionals[END_REF][START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF][START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF][START_REF] Guzina | From imaging to material identification: A generalized concept of topological sensitivity[END_REF][START_REF] Laurain | Topological sensitivity analysis in fluorescence optical tomography[END_REF][START_REF] Masmoudi | The topological asymptotic expansion for the Maxwell equations and some applications[END_REF]. The topological derivative of an objective functional quantifies its perturbation resulting from the virtual creation inside the medium of an object occupying a region B a (z) with prescribed center z and vanishingly small radius a. In the presently-relevant case of a small penetrable obstacle embedded in a three-dimensional acoustic medium, the well known expansion J(a) = J(0) + a 3 T (z) + o(a 3 ) holds [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], where J(a) is the value of the objective functional of interest when the medium hosts a scatterer with support B a (z) and prescribed material parameters. The topological derivative field z → T (z) can then be used as a means for qualitative flaw identification (or, in other contexts, to steer topology optimization algorithms [START_REF] Samet | The topological asymptotic for the Helmholtz equation[END_REF]). Such studies usually involve objective functionals that depend on B a (z) implicitly through the acoustic field u a arising in the perturbed medium, i.e. of the form J(a) = J(u a ). The scattered field v a := u a -u at any given location x = z (where u is the prescribed incident field constituting the applied excitation) is known from many previous studies, e.g. [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF][START_REF] Bendali | Approximation by multipoles of the multiple acoustic scattering by small obstacles in three dimensions and application to the Foldy theory of isotropic scattering[END_REF], to verify v a (x) = O(a 3 ) for three-dimensional configurations.

The topological derivative concept can be extended towards expanding J(a) to higher orders in a, yielding more accurate asymptotic approximations of J. Previous studies following this path include [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of misfit function[END_REF] for sound-hard obstacles in 3D acoustic media, [START_REF] Bonnet | Fast identification of cracks using higher-order topological sensitivity for 2-D potential problems[END_REF] for cracks in 2D conducting media, [START_REF] Ferreira | A new non-iterative reconstruction method for the electrical impedance tomography problem[END_REF] for 2D electrical impedance tomography (EIT). Elastic media are considered in [START_REF] Silva | Higher order topological derivatives in elasticity[END_REF] (expansion of the potential energy in 2D solids with small holes) and [START_REF] Bonnet | Higher order topological derivatives for three-dimensional anisotropic elasticity[END_REF] (general objective functionals for 3D solids containing small inhomogeneities). This kind of higher-order asymptotic approximation may be used as a surrogate of the original functional J for minimization-based identification. This approach, which was in particular shown in [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of misfit function[END_REF][START_REF] Bonnet | Fast identification of cracks using higher-order topological sensitivity for 2-D potential problems[END_REF] to be fast and quantitatively accurate, defines a sampling method insofar as a spatial region of interest may be sampled with a large number of trial flaw sites z, since the minimization of a polynomial approximation of J(a) with respect to a for given z entails very little computational work.

This article, which continues work initiated in [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of misfit function[END_REF][START_REF] Bonnet | Fast identification of cracks using higher-order topological sensitivity for 2-D potential problems[END_REF], addresses the establishment and use of higher-order topological expansions of objective functionals for the case of a small penetrable obstacle embedded in a three-dimensional acoustic medium, under time-harmonic conditions. Such expansions have the form (1) J(a) = J(0) + a 3 T 3 (z) + a 4 T 4 (z) + a 5 T 5 (z) + a 6 T 6 (z) + o(a 6 ) = J 6 (a) + o(a 6 ), where higher-order topological derivatives T 4 (z), T 5 (z), T 6 (z) appear in addition to the previously-known topological derivative T 3 (z) = T (z). The salient features of this work are as follows: (a) The main result is an expansion of the form [START_REF] Ammari | Recovery of small inhomogeneities from the scattering amplitude at a fixed frequency[END_REF] for obstacles of arbitrary shape and a wide class of objective functionals. The case of least-squares (LS) cost functionals, which underpins the most common computational inversion methods, is of primary interest in this work. For this case, the adopted expansion order O(a 6 ) is a natural choice, because a 6 T 6 (z) is the lowest order at which the quadratic term in v a contributes to expansion [START_REF] Ammari | Recovery of small inhomogeneities from the scattering amplitude at a fixed frequency[END_REF]. Evidence from numerical experiments using LS functionals shows that we often have T 6 (z) > 0 (helped by the contribution of the second-order derivative J (u; v a ), which is positive in this case) but T 5 (z) < 0, implying that the O(a 6 ) approximation (1) of J(a) has a minimum for a ∈ R + whereas the O(a 5 ) approximation does not. We both derive the relevant high-order topological expansion (supplementing previously-known results on T 3 (z) by complete expressions for functions T 4 (z), T 5 (z), T 6 (z)) and provide its justification. By contrast, the previous related investigations [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of misfit function[END_REF][START_REF] Bonnet | Fast identification of cracks using higher-order topological sensitivity for 2-D potential problems[END_REF][START_REF] Silva | Higher order topological derivatives in elasticity[END_REF] neither consider penetrable scatterers nor justify the order of approximation achieved by the formally-derived expansions. (b) A computationally light approximate global search procedure (previously introduced in [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of misfit function[END_REF] for sound-hard obstacles, and also used in [START_REF] Ferreira | A new non-iterative reconstruction method for the electrical impedance tomography problem[END_REF] for EIT), where the location and size of the sought object are estimated by minimizing J 6 over a search grid, is presented and demonstrated on numerical experiments for the identification of a penetrable scatterer embedded in a semi-infinite medium. (c) Like in [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of misfit function[END_REF][START_REF] Bonnet | Fast identification of cracks using higher-order topological sensitivity for 2-D potential problems[END_REF], the asymptotic expansion exploits the adjoint solution associated with J, which allows to establish (1) on the basis of (i) the (inner) expansion of u a inside B a and (ii) the leading-order (outer) expansion of u a on the support of the objective function density.

(d) Asymptotic expansions of u a have been the subject of previous studies. Highorder expansions for the Dirichlet problem involving penetrable objects are established in [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF] using coupled boundary integral equation formulations. Matched expansions are used in [START_REF] Ammari | Electromagnetic scattering by small dielectric inhomogeneities[END_REF] (to obtain uniform O(a 3 ) expansions for penetrable objects) and [START_REF] Bendali | Approximation by multipoles of the multiple acoustic scattering by small obstacles in three dimensions and application to the Foldy theory of isotropic scattering[END_REF] (who give arbitrary-order expansions for impenetrable objects). Arbitrary-order expansions are given in [START_REF] Ammari | Polarization and Moment Tensors with Applications to Inverse Problems and Effective Medium Theory[END_REF] for zero-frequency problems in conducting or elastic media. To derive the required solution expansion for the present case (bounded acoustic medium with a penetrable obstacle, Neumann boundary conditions), this work exploits a volume integral equation (VIE) setting [START_REF] Bonnet | A modified volume integral equation for anisotropic elastic or conducting inhomogeneities. Unconditional solvability by Neumann series[END_REF][START_REF] Martin | Acoustic scattering by inhomogeneous obstacles[END_REF], which is natural for modelling many inhomogeneity problems.

The geometrical support of the VIE is B a ; this facilitates the use of coordinate rescaling commonly used in the derivation of asymptotic models and, in combination with the adjoint solution approach, makes the inner expansion in B a play a key role for both the derivation and the justification of the expansion of u a . The adopted VIE setting therefore fits well the present goals, while having broader usefulness as it is transposable to the derivation of similar asymptotic approximations for many other physical contexts, e.g. elasticity [START_REF] Bonnet | Higher order topological derivatives for three-dimensional anisotropic elasticity[END_REF]. The article is organised as follows. Section 2 introduces the acoustic forward and inverse problems, the objective functionals and small-obstacle configurations of interest and the adjoint solution approach. Then, the governing VIE for the forward problem is expanded in Section 3, to obtain the O(a 4 ) inner expansion of u a , the known leading-order outer expansion of u a being recovered afterwards. These results are used in Section 4 to formulate the topological derivatives T 3 , . . . , T 6 and justify the resulting expansion (1) of J(a) for a single penetrable obstacle of arbitrary shape (Theorem 1). The common case of a centrally-symmetric obstacle (Sec. 4.2) is then shown to yield simpler formulas, which become almost explicit for spherical or ellipsoidal shapes. After briefly explaining how Theorem 1 can be generalized to the case of multiple obstacles with fixed locations, and discussing some (e.g. computational) side issues, in Section 5, a simple approximate global search procedure based on the expansion (1) is outlined in Section 6 and demonstrated on numerical experiments in Section 7. Finally, Sections 8 and 9 give the proof of two sub-results (Propositions 1) and 2) upon which Theorem 1 relies.

Acoustic problem and objective functional.

Let Ω ⊂ R 3 be a Lipschitz domain, filled with a homogeneous acoustic medium with mass density ρ 0 and wave velocity c 0 ; this configuration will be called the background (i.e. obstaclefree) medium. All field quantities are assumed to be time-harmonic, with angular frequency ω and the time-harmonic factor e -iωt omitted and implicitly understood throughout. The main exposition is focused on the case where Ω is bounded, its boundary ∂Ω being assumed for definiteness to support prescribed values V D of the normal velocity (straightforward adaptations to the main development allow to consider other types of boundary conditions, see Sec. 5.2 and also Sec. 5.3 for the unbounded case Ω = R 3 ). This excitation gives rise to the background (i.e. incident) acoustic pressure field u, which solves the Neumann boundary-value problem

(2) ∆u + k 2 u = 0 in Ω, ∂ n u = iρ 0 ωV D on ∂Ω,
where k := ω/c 0 is the background wavenumber, n is the normal on ∂Ω outward to Ω, and ∂ n w := ∇w•n is the normal derivative of a field w (the symbol '•' denoting an inner product). It is assumed that ω is not an eigenfrequency of the boundary-value problem [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF] with homogeneous data.

2.1. Forward scattering problem. The incident field u serves to probe Ω for hidden objects. When the medium hosts a penetrable obstacle with support B Ω and homogeneous material parameters (ρ, c), the same applied excitation V D gives rise to a total field u B solving the transmission problem (with the relative material contrasts β, η defined by β := ρ 0 /ρ -1 and η := ρ 0 c 2 0 /ρc 2 -1):

∆u B + k 2 u B = 0 in Ω \ B, (β + 1)∆u B + (η + 1)k 2 u B = 0 in B, ∂ n u B = iρ 0 ωV D on ∂Ω, u B | + = u B | -and ∂ n u B | + = (β + 1)∂ n u B | - on ∂B.
2.2. Inverse problem. We consider the inverse problem of identifying an unknown object ( B, ρ,c) from partial knowledge of the solution ůB of the transmission problem (3) written for ( B, ρ,c), namely values u obs of the acoustic pressure on the measurement surface S obs ⊂ ∂Ω (i.e. u obs = ůB | S obs in the absence of measurement noise). The misfit between the observation u obs and its acoustic prediction u B for a trial obstacle (B, ρ, c) is quantified through an objective functional J (B, β, η) to be minimized. We consider generic objective functionals having the format

(3) J (B, β, η) = J(u B ), with J(w) = ∂Ω ϕ w R (x), w I (x), x dS(x),
where the real-valued density ϕ is a C 2 function of its first two arguments and subscripts 'R', 'I' indicate the real and imaginary parts of a complex number, i.e. w R = Re(w), w I = Im(w). We also denote by ∂ R ϕ and ∂ I ϕ the partial derivatives of ϕ with respect to its first and second argument, respectively. The corresponding second-order derivatives of ϕ, similarly denoted (e.g. ∂ RI ϕ), are additionally assumed to be C 0,γ functions of their first two arguments for some γ > 0. Data u obs recorded on S obs ⊂ ∂Ω is easily accommodated in (3); for instance, for a leastsquares misfit functional that is customarily applied to inverse problems, we have (with χ A denoting the characteristic function of a set A)

(4) ϕ(w R , w I , •) = 1 2 w -u obs 2 χ S obs = 1 2 (w R -u obs R ) 2 + (w I -u obs I ) 2 χ S obs , Remark 1.
The chosen nature of the data u obs is for definiteness, other types of data requiring only straightforward adjustments to this setting and the forthcoming main development, see Sec. 5.2.

2.3.

Identification by objective functional expansion. Letting B ⊂ R 3 denote a fixed bounded open set containing the origin, with volume |B|, we consider a penetrable obstacle of support B a (z) = z + aB and size a > 0, endowed with relative material contrasts (β, η) and centered at a given point z ∈ Ω. Without loss of generality, z is assumed to be the centroid of B a . We then set B = B a (z) in the transmission problem (3) and objective functional [START_REF] Ammari | Polarization and Moment Tensors with Applications to Inverse Problems and Effective Medium Theory[END_REF]. Denoting by u a (•; z) the total field solving problem (3), we define the cost function J(a; z) in terms of J by ( 5)

J(a; z) = J(u a ) = ∂Ω ϕ u a R (x), u a I (x), x dV (x)
(for notational convenience, explicit references to z will sometimes be omitted in the sequel, e.g. by writing J(a) or u a (x) instead of J(a; z) or u a (x; z)). The inverse scattering problem is then recast in terms of finding a, z that minimize an asymptotic expansion of J(a; z) in powers of a, i.e. a polynomial in a with coefficients depending on z (and also on B, β, η) that approximates J(a; z).

2.4.

Cost functional expansion using adjoint solution. The desired asymptotic expansion of J(a) is sought on the basis of an expansion of J(u a ) to second order about u a = u, i.e.:

(6) J(a) = J(u) + J (u; v a ) + 1 2 J (u; v a ) + R(u; v a ), where v a := u a -u is the scattered field caused by the small obstacle, and

J (u; w) = p=R,I ∂Ω ∂ p ϕ u R , u I , • w p dS, J (u; w) = p,q=R,I ∂Ω ∂ pq ϕ u R , u I , • w p w q dS, (7) R(u; w) = t 0 (1 -t)J (u + tw; w) dt -1 2 J (u; w).
To evaluate J (u; v a ), it is convenient to introduce the adjoint field û, defined as the solution of the problem

(8) ∆û + k 2 û = 0 in Ω, ∂ n û = -∂ R ϕ -i∂ I ϕ (u R , u I , •) on ∂Ω,
whose boundary data is chosen such that J (u; w) = -Re ∂ n û, w ∂Ω . Then, J (u; v a ) can be expressed in terms of the restrictions of u and û on the small obstacle support B a . For convenience in this work, and in particular to express J (u; v a ) in concise form, we introduce the bilinear forms

(9) u, w R X := X R ∇u•∇w dV, u, w E X := X E uw dV,
defined on H 1 (X) × H 1 (X), for a domain X ⊆ Ω and coefficients R, E ∈ L ∞ (X), and additionally set u, w X := u, w R=1 X and u, w X := u, w E=1 X , with similar notation used for surface integrals if appropriate. Using these notations, we have: Lemma 1. Let u a and û be the total and adjoint fields, respectively solving problems (3) with B = B a and (8). The first-order derivative J (u; v a ) in (6) is then given by

J (u; v a ) = Re u a , û β Ba -k 2 u a , û η Ba .
Proof. Writing in weak form the adjoint problem (8) (with test function v a ) and the forward problems (2) and (3) (with test function û), we have the identities

û, v a Ω -k 2 û, v a Ω = -∂ R ϕ -i∂ I ϕ, v a ∂Ω , (a) u, û Ω -k 2 u, û Ω = iρ 0 ω V D , û ∂Ω , (b) u a , û Ω -k 2 u a , û Ω + u a , û β Ba -k 2 u a , û η Ba = iρ 0 ω V D , û ∂Ω (c).
The lemma results from the combination -(a)-(b)+(c) and definition (7) of J .

Remark 2. For the least-squares objective functional defined by (4), we have R(u; v a ) = 0 and

J (u; w) = S obs
Re (u -u obs )w dS, J (u; w) = S obs |w| 2 dS.

2.5. Summary of previous results on topological derivative. The leading contribution to J(u a ) as a → 0 is given by that of J (u, v a ), which 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] as

(10) J (u, v a ) = a 3 T 3 (z) + o(a 3 )
in terms of the topological derivative T 3 (z), given (for the present physical context by

T 3 (z) = -Re ∇û•A(β) 11 •∇u -|B| k 2 ηûu (z)
, where the second-order polarization tensor A 11 has been established for any obstacle shape B in e.g. [START_REF] Ammari | Polarization and Moment Tensors with Applications to Inverse Problems and Effective Medium Theory[END_REF][START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF][START_REF] Ammari | Recovery of small inhomogeneities from the scattering amplitude at a fixed frequency[END_REF] (see e.g. equation 56 when B is an ellipsoid). Moreover, the leading asymptotic behaviour of the total field inside B a and on the external surface is respectively characterized by

u a (x; z) = u(z) + aU 1 (x -z)/a + o(a)
x ∈ B a , (11)

u a (x; z) = u(x) + a 3 W (x; z) + o(a 3 ) x ∈ ∂Ω, (12)
where U 1 , W are known functions that depend on B and β, η (see equations (39b) and (50)).

2.6. Deriving the expansion of J(a): preliminary considerations and methodology. We seek an expansion of J(u a ) that includes the leading contribution as a → 0 of J (u, v a ). In view of ( 6) and ( 12), the lowest expansion order meeting this goal is O(a 6 ). Moreover, the known outer solution expansion [START_REF] Bonnet | Higher order topological derivatives for three-dimensional anisotropic elasticity[END_REF] suffices for determining the (leading) O(a 6 ) expansion of J (u; v a ) through [START_REF] Bellis | Acoustic inverse scattering using topological derivative of far-field measurements-based L 2 cost functionals[END_REF], so that the main task at hand is to find the O(a 6 ) expansion of J (u; v a ). Lemma 1 shows that this requires expanding u a in the vanishing obstacle B a . To this aim, we recast the forward problem (3) as a volume integral equation (VIE) whose geometrical support is B a , making the restriction of u a to B a its main unknown. Then, following the usual approach for solution asymptotics involving a small length parameter, the position vector x in B a is scaled according to x = z + ax, and a rescaled main unknown U a such that U a (x) := u a (z+ax) is introduced. An asymptotic expansion of U a is then sought in the form U a = U 0 + aU 1 + . . ., consistently with [START_REF] Bonnet | A modified volume integral equation for anisotropic elastic or conducting inhomogeneities. Unconditional solvability by Neumann series[END_REF]. On expressing J (u; v a ) as given by Lemma 1 in terms of U a and integrals over B, noting that the coordinate stretching [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] implies ∇u a (x) = a -1 ∇U a (x), we obtain

(13) J (u; v a ) = a 3 Re û, U a •-z a β Ba -k 2 û, U a •-z a η
Ba . Expanding J (u; v a ) to order O(a 6 ) therefore entails finding the O(a 4 ) expansion of U a in B. The O(1) term of the expansion of U a will be found to be constant in B (as physical intuition suggests), so that ∇U a = O(a) and ( 13) is consistent with the expected leading-order expansion [START_REF] Bonnet | Fast identification of cracks using higher-order topological sensitivity for 2-D potential problems[END_REF].

We carry out the above-outlined plan in the next parts: the needed elements of the expansion of u a are derived and justified next in Section 3; then, the O(a 6 ) expansion of J(a) is established in Section 4.

3.

Inner and outer expansions of the total field. Focusing on u a in B a as the primary quantity, we employ a volume integral equation (VIE) formulation (Sec. 3.1), express it using scaled coordinates (Sec. 3.2) and derive its leading asymptotic form (Sec. 3.3). This reveals the essential role played by a family of zero-frequency (Laplace) free-space transmission problems (FSTPs), surveyed in Sec. 3.4. Finally, the inner solution expansion is carried out in Section 3.5, and the leading outer solution expansion recalled in Sec. 3.6.

3.1.

Volume integral equation formulation for the forward problem. We begin by stating the relevant VIE formulation. Let G be the Green's function associated with the domain Ω with a Neumann boundary condition, solving the BVP [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF] (a) ∆G(•,

x) + k 2 G(•, x) = δ x in Ω, (b) ∂ n G(•, x) = 0 on ∂Ω
(with δ x denoting the Dirac distribution supported at x). Then, the total field u a inside B a associated with the solution of (3) with B = B a is governed by the volume integral equation (VIE)

(15) (I -L a )u a = u in B a ∪ (Ω \ B a )
in which I is the identity operator and the integral operator L a is defined by

(16) L a w(x) = k 2 Ba η G(y, x)w(y) dV (y) - Ba β ∇ 1 G(y, x)•∇w(y) dV (y)
where ∇ k G denotes the gradient with respect to its k-th argument of a function G having two (or more) arguments; likewise, ∇ k G will denote the second-order gradient of G with respect to its k-th and -th arguments. Known mapping properties of integral operators treated as pseudodifferential operators [22, Thm. 6.1.12] imply that L a is well defined as a H 1 (B a ) → H 1 (B a ) operator. A concise derivation of the VIE ( 15) is given in Appendix A; otherwise see [START_REF] Bonnet | A modified volume integral equation for anisotropic elastic or conducting inhomogeneities. Unconditional solvability by Neumann series[END_REF] and references therein. In the case Ω = R 3 , for which G = G ∞ k (see eq. ( 23)), the VIE ( 15) is identical (up to notational adjustment) to that established in [START_REF] Martin | Acoustic scattering by inhomogeneous obstacles[END_REF] and applied to the obstacle configuration (B a , β, η). 15) is now recast in rescaled form (i.e. using integral operators defined on the fixed compact domain B) by effecting the coordinate change [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] x = z + ax

Rescaled form of VIE formulation. The integral equation (

(x ∈ B a , x ∈ B),
noting that this process also transforms the differential volume element according to dV x = a 3 d Vx (x ∈ B a , x ∈ B). Towards this reformulation, we introduce the mapping P a : H 1 (B a ) → H 1 (B) that applies the coordinate scaling [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] to functions, and its inverse P -1 a : (18)

P a v(x) := v(z + ax) , P -1 a V (x) := V (x -z)/a .
The coordinate stretching exerted by [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] then implies that gradients transform under P a according to

(19) P a [∇v](x) = a -1 ∇(P a v)(x)
Applying the scaling (17) to both arguments x, y in the VIE [START_REF] Cedio-Fengya | Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction[END_REF], we obtain the governing VIE (20

) (I -L a )U a = P a u in B, with L a = P a L a P -1 a . for the rescaled unknown (21) U a (x) := P a u a (x) = u a (z + ax) x ∈ B,
where the scaled version

L a := P a L a P -1 a : H 1 (B) → H 1 (B) of the integral operator L a is given by L a U (x) = k 2 a 3 B η G z + aȳ, z + ax U (ȳ) dV (ȳ) (22) -a 2 B β ∇ 1 G z + aȳ, z + ax •∇U (ȳ) dV (ȳ)
3.3. Asymptotic behavior of the integral operator. We now seek the asymptotic form of the operator L a as a → 0. This entails expanding the Green's function G. To this aim, the decomposition

(23) G(y, x) = G ∞ 0 (y -x) + G ∞ k (y -x) -G ∞ 0 (y -x) + G c (y, x) of G is used, with G ∞
k and G ∞ 0 , the singular fundamental solutions of the Helmholtz and Laplace equations, given by

G ∞ k (r) = e ik|r| 4π|r| , G ∞ 0 (r) = 1 4π|r|
.

The complementary acoustic Green's function G c (•, x) is seen (using superposition) to solve a BVP of the form (2) with boundary condition

∂ n G c = -∂ n G ∞ k (• -x) on ∂Ω; therefore, G c (•, x) ∈ C ∞ (Ω) and is in particular bounded at x. Since both G ∞ k -G ∞ 0 and ∇G ∞ k -∇G ∞ 0 are bounded at r = 0, the singular parts of G and ∇ 1 G are entirely contained in G ∞
0 and ∇G ∞ 0 , respectively. Moreover, G ∞ 0 is positively homogeneous with degree -1, so that the coordinate scaling [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] implies

G ∞ 0 (y -x) = a -1 G ∞ 0 (ȳ -x), ∇G ∞ 0 (y -x) = -a -2 ∇G ∞ 0 (ȳ -x).
The above properties, combined with [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF], allow to recast the scaled integral operator L a given by ( 22) as [START_REF] Laurain | Topological sensitivity analysis in fluorescence optical tomography[END_REF] L a = H + a 2 L c a , where the leading operator H and the complementary operator L c a are defined by

HU (x) = B β ∇G ∞ 0 (x -•)•∇U dV L c a U (x) = k 2 a B η G z + aȳ, z + ax U (ȳ) dV (ȳ) - B β ∇ 1 G c z + aȳ, z + ax •∇U (ȳ) dV (ȳ) - B β ∇ G ∞ k -G ∞ 0 (a(ȳ -x) •∇U (ȳ) dV (ȳ), (25)
noting that H does not depend on a. The VIE [START_REF] Guzina | Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics[END_REF] for U a then takes the form

(26) (I -H -a 2 L c a )U a = P a u in B,
where I -H is the leading term of L a as a → 0. In fact, the scaled operators I -H and L c a and their original (unscaled) counterparts P -1 a I -H P a and L c a := a 2 P -1 a L c a P a have the following properties:

Proposition 1. Let β such that -1 < β < ∞. 1. The operator I -H : H 1 (B) → H 1 (B) is (a) well defined and bounded, (b) boundedly invertible. 2. Let a 0 such that B a Ω for any a ≤ a 0 . Then: (i) The operator A a := P -1 a I -H P a : H 1 (B a ) → H 1 (B a
) is well-defined, bounded and invertible. Moreover, there exists C 0 > 0 such that A -1 a ≤ C 0 for any a < a 0 . (ii) There exists C c > 0 such that, for any a < a 0 , the operator L c a :

H 1 (B) → H 1 (B) satisfies L c a ≤ C c and the operator L c a := a 2 P -1 a L c a P a : H 1 (B a ) → H 1 (B a ) satisfies L c
a ≤ aC c . Regarding item 1, (a) stems from Thm. 6.1.12 in [START_REF] Hsiao | Boundary Integral Equations[END_REF] on pseudodifferential operators, and (b) is [START_REF] Bonnet | A modified volume integral equation for anisotropic elastic or conducting inhomogeneities. Unconditional solvability by Neumann series[END_REF]Thm. 3.1]. The proof of item 2 is given in Section 8. Remark 3. We note that the non-leading parts of L a and L a are found to have respective orders O(a) and O(a 2 ) (the leading order being O(1) for both the original and the rescaled VIEs); this order disparity is due to the effect of the coordinate stretch [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] on density gradients, see [START_REF] Gumerov | Fast multipole methods for the Helmholtz equation in three dimensions[END_REF].

The rescaled version [START_REF] Masmoudi | The topological asymptotic expansion for the Maxwell equations and some applications[END_REF] of the VIE [START_REF] Cedio-Fengya | Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction[END_REF] provides the basis for formally deriving the inner expansion of u a . Since I -H is the leading part of the integral operator I-L a (see Prop. 1), VIEs of the form (I-H)u B = u, which in fact govern (frequencyindependent) free-space transmission problems (FSTPs), will play a key role, as coefficients of the expansion of u a will be found to satisfy such VIEs. Accordingly, some definitions and results for FSTPs are now gathered for later use.

3.4.

Free-space transmission problems. Integral equations of the form

(27) (I -H)u B = u in R 3 ,
are known to govern a Laplace free-space transmission problem (FSTP) for an inhomogeneity B with (e.g. conductivity) relative contrast β embedded in an infinite medium, defined as follows (see e.g. [START_REF] Bonnet | A modified volume integral equation for anisotropic elastic or conducting inhomogeneities. Unconditional solvability by Neumann series[END_REF]): given a background field u ∈ H 1 loc (R 3 ), find the total field u B such that [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] div

(1 + βχ B )∇u B = div (∇u) in R 3 , u B (y) -u(y) = O(|y| -2 ) (|y| → ∞).
Solutions to such FSTPs satisfy the following reciprocity identity (see proof in Appendix D.1):

Lemma 2. Let u B , u B solve (28) for respective background fields u, u . We have u B , u

β B = u B , u β B .
Polynomial background field. The case where the background field u is polynomial will play an important role. Any polynomial of degree n may be set in the form [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF] u

(x) = E 0 + n p=1 π p [E p ](x) with π p [E p ](x) := E p • x ⊗p ,
where E 0 ∈ C is a scalar and E p are constant complex tensors of order p that are assumed without loss of generality to be invariant under all index permutations. In [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF] and hereinafter, the notation A • B indicates, for any tensors A, B of respective orders p, q, their m-fold inner product, with m := min(p, q) (i.e. the inner product is taken over as many indices as possible), e.g. (A • B) ijk = A ijk mn B mn (using Einstein's summation convention); moreover, the notation x ⊗p is a shorthand for the tensor product x ⊗ . . . ⊗ x of order p (e.g.

x ⊗3 = x ⊗ x ⊗ x, i.e. (x ⊗3 ) ijk = x i x j x k ).
The above definitions, and in particular this convention, imply that ∇π 1) . We then denote by U p B [E p ] the solution to the FSTP [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] with u = π p [E p ], which satisfies the well-posed integral equation [START_REF] Samet | The topological asymptotic for the Helmholtz equation[END_REF] (

p [E p ](x) = pE • x ⊗(p-
I -H)U p B = π p [E p ] in B. Since H[E 0 ] = 0, U 0 B is a constant function: U 0 B = E 0 . Moreover, when B is an ellip- soid, it is well known that U p B [E p ] is polynomial with degree p in B (Appendix B
). We will need cases p = 1, 2, for which

(31) ∇U 1 B (y) = [ I + βS 1 ] -1 •E 1 , ∇U 2 B (y) = [ I + βS 2 ] -1 : E 2 •y y ∈ B
(with I and I denoting the second-order identity tensor and the fourth-order identity for symmetric second-order tensors, respectively), where S 1 and S 2 are constant tensors, respectively of order 2 and 4 (see Appendix B). Analytical expressions of S 1 and S 2 are available; in particular, if B is a ball, we have (Appendix C) (32) (a) S 1 = 1 3 I, (b) S 2 = 2 5 2I + I ⊗ I . Polarization tensors (PTs). They appear in many asymptotic expansions involving small inhomogeneities, see e.g. [START_REF] Ammari | Polarization and Moment Tensors with Applications to Inverse Problems and Effective Medium Theory[END_REF][START_REF] Cedio-Fengya | Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction[END_REF]. For integers p, q ≥ 1, let the PT A pq be the constant tensor of order p + q such that (33)

U p B [E p ] , π q [E q ] β B = E p • A pq •
E q for any constant tensors E p and E q of respective orders p, q (equation (33) uniquely defines A pq under the additional condition that A pq be invariant under any permutation of either its first p indices or its last q indices, reflecting the assumptions previously made on E p , E q ). The practical evaluation of tensors A pq is made easier by applying Lemma 2 to (u

B , u B ) = U p B [E p ], U q B [E q ] and (u, u ) = π p [E p ], π q [E q ] , yielding (34) E p • A pq • E q = E q • A qp • E p .
Consequently, knowing the FSTP solution U m B [E m ] with m := min(p, q) is sufficient for evaluating A pq .

Remark 4. The present PT A pq can be shown to define the same bilinear form E p • A pq • E q as its counterpart m ij of [START_REF] Ammari | Polarization and Moment Tensors with Applications to Inverse Problems and Effective Medium Theory[END_REF] (defined in terms of layer potential representations for U p B [E p ], and where i, j are multi-indices of respective lengths |i| = p, |j| = q) whenever the polynomial π p [E p ] is harmonic.

Inner solution expansion.

3.5.1. Formal derivation. Recalling the rescaling [START_REF] Guzina | From imaging to material identification: A generalized concept of topological sensitivity[END_REF], the desired inner expansion of u a in B a is assumed to have the form

u a (x) = P -1 a [U (4) a ](x) + δ a (x) (x ∈ B a ), (35)
with U (4) a := U 0 + aU 1 + a 2 U 2 + a 3 U 3 + a 4 U 4 in terms of functions U 0 , . . . U 4 defined in B, and where the remainder δ a is small in a sense made precise later (see Prop. 2). We now seek governing equations for U 0 , . . . U 4 by substituting the ansatz (35) into the asymptotic expansion form [START_REF] Masmoudi | The topological asymptotic expansion for the Maxwell equations and some applications[END_REF] of the VIE, i.e. by writing that the equation

(36) I -H -a 2 L c a U (4) a = P a u in B. is verified at orders O(1)
to O(a 4 ) included. The desired equations for U 0 , . . . U 4 then result from term-by-term identification of the expansion of the above equation in powers of a, starting from the lowest O(1) order. The right-hand side is expressed by its Taylor expansion (37)

P a [u](x) = u z + 4 p=1 a p p! π p [g p z ](x) + o(a 4 ) =: T 4 [u](x) + o(a 4 )
(having introduced the convenient shorthand notations u z := u(z) and g k z := ∇ k u(z) for the background field and its derivatives at z, and with monomials π p as defined by ( 29

I -H U 1 = π 1 [g z ],
with g z := g 1 z . Having the form [START_REF] Samet | The topological asymptotic for the Helmholtz equation[END_REF], respectively with p = 0 and p = 1, their solutions are (see Sec. 3.4)

(39) (a) U 0 = u z , (b) U 1 = U 1 B [g z ] in B.
We note in particular that ∇U 0 = 0, which will bring some simplification in forthcoming expansions and shows that the leading contribution of ( 13) is, as expected (see Sec. 2.5), of order O(a 3 ).

(ii) Expansion of L c a U (4)
a . Using Taylor expansions of a → G c z + aȳ, z + ax (which is smooth in a neighbourhood of a = 0) and a → e ika|ȳ-x| in [START_REF] Martin | Acoustic scattering by inhomogeneous obstacles[END_REF], the expansion of

L c a U for given U ∈ H 1 (B) is L c a U = L c 0 U + aL c 1 U + a 2 L c 2 U + o(a 2 ), (40) L c 0 U (x) := k 2 (G -N )U - B β∇ 1 G c z •∇U dV L c 1 U (x) := k 2 ik 4π + G c z B η U dV - B β ∇ 11 G c z • ȳ + ∇ 12 G c z • x + ik 3 12π r •∇U dV L c 2 U (x) := -k 2 MU + B η ∇ 2 G c z • x + ∇ 1 G c z • ȳ U dV
(with r := xȳ) in terms of the integral operators G, N and M defined by ( 41)

GU (x) = B η 4π|r| U dV, N U (x) = B β 8π ∇|r|•∇U dV, MU (x) = B η 8π |r|U dV.
We note that all terms in formulas (40) except those involving operators M, N , G yield constant or linear functions of x for any U . The above definitions then result in the following expansion of L c a U

a :

L c a U (4) a = L c 0 U 0 + a(L c 0 U 1 + L c 1 U 0 ) + a 2 (L c 0 U 2 + L c 1 U 1 + L c 2 U 0 ) + o(a 2 ) = k 2 u z G1 + a k 2 (G -N )U 1 + D 3 z + a 2 k 2 (G -N )U 2 -k 4 u z M1 + D 4 z + E z • x + o(a 2
), wherein the scalar constants D 3 z , D 4 z ∈ C and the constant vector E z = C 3 can be evaluated by using (39) and expressing relevant integrals on B in terms of the polarization tensors defined by (33) and (46a) whenever possible. Introducing for convenience the shorthand notations

G c z := G c (z, z), ∇ k G c z := ∇ k G c (z, z) and ∇ k G c z := ∇ k G c (z, z), this evaluation yields (42) D 3 z = k 2 u z η|B| ik 4π + G c z -g z •A 11 •∇ 1 G c z , E z = k 2 u z η|B| ∇ 2 G c z -g z •A 11 • ∇ 12 G c z + ik 3 12π I ,
whereas the precise value of D 4 z will prove irrelevant. (iii) Determination of U 2 , U 3 , U 4 . The O(a 2 ), O(a 3 ) and O(a 4 ) equations arising from (36) are

(a) I -H U 2 = 1 2 π 2 [g 2 z ] + L c 0 U 0 , (b) I -H U 3 = 1 6 π 3 [g 3 z ] + L c 0 U 1 + L c 1 U 0 , ( 43
) (c) I -H U 4 = 1 24 π 4 [g 4 z ] + L c 0 U 2 + L c 1 U 1 + L c 2 U 0 . Their solutions are (44) (a) U 2 = 1 2 U 2 B [g 2 z ]+X B , (b) U 3 = 1 6 U 3 B [g 3 z ]+Y B , (c) U 4 = 1 24 U 4 B [g 4 z ]+Z B with the auxiliary H 1 (B) functions X B , Y B , Z B defined as the solutions of the FSTPs (a) I -H X B = k 2 u z G1 (b) I -H Y B = k 2 (G -N )U 1 + D 3 z (45) (c) I -H Z B = k 2 (G -N )U 2 -k 4 u z M1 + D 4
z + E z . As usual in asymptotic methods, the U are defined sequentially and depend on lower-order solutions.

The following lemma, whose proof is given in Appendix D, gathers identities verified by X B , Y B , Z B . They involve additional constant tensors B pq and Q pq defined, in terms of the FSTP solutions U p B , by

E p • B pq • E q = U p B [E p ] , π q [E q ] η B , (46a)
E p • Q pq • E q = U p B [E p ] , U q B [E q ] η B + U p B [E p ] , N U q B [E q ] β B . ( 46b
)
where N is the operator defined in (41). We note that B 01 = 0 (due to U 0 B being a constant function and z being the centroid of B by assumption). Lemma 3. The functions X B , Y B , Z B solving the FSTPs (45) satisfy the identities

X B , π q [E q ] β B = k 2 u z B 0q • E q -E q • B q0 q ≥ 1 Y B , π q [E q ] β B = k 2 g z •(B 1q -Q 1q )•E q q ≥ 1 Z B , π 1 [ĝ z ] β B = k 2 1 2 g 2 z •(B 21 -Q 21 )•ĝ z + X B , π 1 [ĝ z ] η B -u z 1, Y B η B + u z η|B| D 3 z + E z •A 11 •g z
with the tensors B pq , Q pq as defined by (46a,b) and where D 3 z and Y B (in the last identity) are defined by (42) and problem (45b) with u z , g z replaced with ûz , ĝz .

Resulting solution expansion and its justification.

To assess the approximation order achieved by (35), we need to estimate the H 1 (B a ) norm of δ a := u a -P -1 a U

a , and in particular to show that its order in a is higher than that of the highest-order term P -1 a U 4 . We have a 4 U 4 H 1 (B) = O(a 4 ). Then, we observe that the rescaling mapping P a (see [START_REF] Ferreira | A new non-iterative reconstruction method for the electrical impedance tomography problem[END_REF]) implies

(47) (a) v L 2 (Ba) = a 3/2 P a v L 2 (B) , (b) ∇v L 2 (Ba) = a 1/2 ∇(P a v) L 2 (B)
for any v ∈ H 1 (B a ). Consequently,

u 4 := P -1 a [a 4 U 4 ] is such that u 4 L 2 (Ba) = a 11/2 U 4 L 2 (B) and ∇u 4 L 2 (Ba) = a 9/2 ∇U 4 L 2 (B)
, and hence can be estimated as u 4 H 1 (Ba) = O(a 9/2 ). To justify the approximation (35), we therefore need to prove that δ a H 1 (Ba) = o(a 9/2 ). The following proposition fulfills this objective: Proposition 2. Let β be such that -1 < β < ∞. There exists a 1 > 0 and a constant C > 0 independent of a such that the expansion error δ a := u a -P -1 a U

(4) a of the inner expansion (35) verifies the estimate

(48) δ a H 1 (Ba) ≤ Ca 11/2
for all a < a 1 .

Proof. Inserting u a given by (35) in the VIE [START_REF] Cedio-Fengya | Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction[END_REF], the inner expansion error δ a satisfies the VIE (

I -L a )δ a = γ a with γ a = u -(I -L a ) P -1 a U (4) a
. The proof of estimate (48) then consists in (i) showing that the inverse of the operator (I -L a ) : H 1 (B a ) → H 1 (B a ) exists and is bounded independently of a for any small enough a and (ii) estimating γ a H 1 (Ba) . This proof is given in Sec. 9; it rests upon the bounded invertibility of I -H, which is ensured by the assumption on β (see Prop. 1).

3.6. Outer expansion. We now turn to the expansion of v a (x) for x / ∈ B a . In this case, v a (x) is given by ( 15) used as an integral representation:

(49) v a (x) = L a [u a ](x) x ∈ Ω \ B a .
Since x ∈ B a , y → G(y, x) is smooth for y ∈ B a , and one therefore has

∇ 1 G(y, x) = ∇ 1 G(z, x)+O(a)
, while the inner expansion (35) truncated to its leading order gives

u a (y) = u(z) + O(a) and ∇u a (y) = g z +P -1 [∇U 1 ](ȳ) + O(a).
Using this together with the rescaling ( 17) in (49) yields the outer expansion

v a (x) = a 3 W (x; z) + O(a 4 ) x ∈ Ω \ B a , ( 50
) with W (x; z) := -∇ 1 G(z, x)•A 11 •g z + k 2 η|B|G(z, x)u z
which holds pointwise and is well-known, e.g. [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF][START_REF] Ammari | Electromagnetic scattering by small dielectric inhomogeneities[END_REF]. Expansion (50) holds pointwise up to the boundary if ∂Ω is C 1,1 (which ensures continuity up to ∂Ω of x → G(z, x) = G(x, z), by e.g. [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Equations[END_REF]Thm. 4.18]). Such extra smoothness for ∂Ω is avoided by the following version of the outer expansion:

Misfit function expansion.

Exploiting the results established and facts gathered in Sec. 3, we now formulate and justify the O(a 6 ) expansion of J(a). The expansion is first given in its most general form, valid for a small obstacle of arbitrary shape, in Theorem 1; this is the main result of this article. Then, the expansion is shown to take a significantly simpler form if restricted to the sub-class of centrally-symmetric inclusions (Sec. 4.2), which includes the important special cases of ellipsoidal and spherical obstacles for which explicit forms of the expansion coefficients can be given (Sec. 4.3).

Expansion in the general case.

Theorem 1. For a single obstacle, characterized by its geometrical support B a := z+aB and relative material parameters β, η and embedded in the three-dimensional homogeneous reference medium Ω, the O(a 6 ) expansion of any objective function J(a) of format (5) with a density ϕ(w R , w I , x) that is twice differentiable in its first two arguments, with C 0,γ second-order derivatives for some γ > 0, is

J(a) = J(0) + a 3 T 3 (z) + a 4 T 4 (z) + a 5 T 5 (z) + a 6 T 6 (z) + o(a 6 )
with the topological derivatives T 3 , . . . T 6 given by

T 3 (z) = Re g z •A 11 •ĝ z -k 2 η|B|u z ûz (51a) T 4 (z) = Re 1 2 A 12 • (g z ⊗ ĝ2 z + ĝz ⊗ g 2 z ) -k 2 B 10 •(u z ĝz + ûz g z ) (51b) T 5 (z) = Re 1 6 A 13 • (g z ⊗ ĝ3 z + ĝz ⊗ g 3 z ) + 1 4 g 2 z : A 22 : ĝ2 z (51c) -1 2 k 2 B 20 : (u z ĝ2 z + ûz g 2 z ) -k 2 g z •Q 11 •ĝ z -k 2 ûz 1, X B η B T 6 (z) = Re 1 2 J (u; W ) + 1 24 A 14 • (g z ⊗ ĝ4 z + ĝz ⊗ g 4 z ) (51d) + 1 12 A 23 • (g 2 z ⊗ ĝ3 z + ĝ2 z ⊗ g 3 z ) -1 6 k 2 B 30 • (u z ĝ3 z + ûz g 3 z ) -1 2 k 2 Q 12 • (g z ⊗ ĝ2 z + ĝz ⊗ g 2 z ) + E z •A 11 •ĝ z -k 2 ûz 1, Y B η B -k 2 1, Y B η B + k 2 |B|η D 3 z ûz .
The tensors A pq , B pq and Q pq are respectively defined by (33), (46a) and (46b) in terms of solutions to free-space transmission problems (FSTPs) with polynomial background field (see Sec. 3.4), the function W is given by (50), the functions X B , Y B solve the FSTPs (45a,b), the scalar and vector constants D 3 z , E z are as in Lemma 3, and the shorthand notations g k z , ĝk z for derivatives at z of u, û are as in (37). Proof. The proof results from expanding in powers of a each term of expansion (6) of J(u a ).

(a) First term of (6). Recalling Lemma 1 and eqn. [START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF], and effecting the coordinate change [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] in the relevant integrals, J (u; v a ) is given by

J (u; v a ) = Re a P a û, U (4) a β B -k 2 a 3 P a û, U (4) a η B (52a) +Re û, δ a β Ba -k 2 û, δ a η Ba
where, inserting (39) and ( 44) into (35), the O(a 4 ) approximation . We first evaluate the term P a û, U (4) a β B . Using (52b) and a Taylor expansion in a of P a û (i.e. (37) with u replaced with û), introducing the polarization tensors A pq and invoking the symmetry property (34) wherever possible, we obtain

U (4) a of U a is given by (52b) U (4) a = u z + aU 1 B + a 2 1 2 U 2 B + X B + a 3 1 6 U 3 B + Y B + a 4 1 24 U 4 B + Z B and that of ∇U
a P a û, U (4) a β B = a 3 g z •A 11 •ĝ z + a 4 1 2 A 12 • (g z ⊗ ĝ2 z + ĝz ⊗ g 2 z ) + 2 X B , π 1 β B + a 5 1 6 A 13 • (g z ⊗ ĝ3 z + ĝz ⊗ g 3 z ) + 1 4 g 2 z : A 22 : ĝ2 z + 1 2 X B , π 2 β B + Y B , π 1 β B + a 6 1 24 A 14 • (g z ⊗ ĝ4 z + ĝz ⊗ g 4 z ) + 1 12 A 23 • (g 2 z ⊗ ĝ3 z + ĝ2 z ⊗ g 3 z ) + 1 6 X B , π 3 β B + 1 2 Y B , π 2 β B + Z B , π 1 β B + O(a 7 ),
having introduced the shorthand notation π m := π[ĝ m z ]. Then, all terms involving auxiliary FSTP solutions X B , Y B , Z B are transformed with the help of Lemma 3, to obtain

P a û, U (4) a β B = a 3 g z •A 11 •ĝ z + a 4 1 2 A 12 • (g z ⊗ ĝ2 z + ĝz ⊗ g 2 z ) -2k 2 ĝz •B 10 u z + a 5 1 6 A 13 • (g z ⊗ ĝ3 z + ĝz ⊗ g 3 z ) + 1 4 g 2 z : A 22 : ĝ2 z + 1 2 k 2 u z B 02 : ĝ2 z -ĝ2 z : B 20 + k 2 g z •(B 11 -Q 11 )•ĝ z + a 6 1 24 A 14 • (g z ⊗ ĝ4 z + ĝz ⊗ g 4 z ) + 1 12 A 23 • (g 2 z ⊗ ĝ3 z + ĝ2 z ⊗ g 3 z ) + E z •A 11 •ĝ z + 1 6 k 2 u z B 03 • ĝ3 z -ĝ3 z • B 30 -1 2 k 2 Q 12 • (g z ⊗ ĝ2 z + ĝz ⊗ g 2 z ) + 1 2 k 2 g z •B 12 : ĝ2 z + 1 2 k 2 g 2 z : B 21 •ĝ z + k 2 X B , π 1 η B -k 2 u z 1, Y B η B + k 2 u z η|B| D 3 z + O(a 7 ) (52c)
Similarly evaluating the term P a û, U (4) a η B by using the inner expansion (52b), we find

P a û, U (4) a η B = a 3 |B|ηu z ûz + a 4 g z •B 10 ûz + a 5 1 2 g 2 z : B 20 ûz + g z •B 11 •ĝ z + 1 2 u z B 02 : ĝ2 z + ûz 1, X B η B + a 6 1 6 g 3 z • B 30 ûz + 1 2 g 2 z : B 21 •ĝ z + 1 2 g z •B 12 : ĝ2 z + 1 6 u z B 03 • ĝ3 z + ûz 1, Y B η B + X B , π 1 η B + O(a 7 ) (52d)
Moreover, the second term in the right-hand side of (52a) admits the following estimate by virtue of Proposition 2 and the Cauchy-Schwarz inequality:

û, δ a β Ba -k 2 û, δ a η Ba (52e) ≤ |β| ∇û L 2 (Ba) ∇δ a L 2 (Ba) + |η| û L 2 (Ba) δ a L 2 (Ba) ≤ |β| + |η| û H 1 (Ba) δ a H 1 (Ba) ≤ Ca 3/2 a 11/2 = Ca 7 .
Combining (52c), (52d) and (52e), we finally find [START_REF] Bellis | A FEM-based topological sensitivity approach for fast qualitative identification of buried cavities from elastodynamic overdetermined boundary data[END_REF]. It admits, as a direct consequence of the outer expansion (50), the expansion

J (u; v a ) = a 3 g z •A 11 •ĝ z -k 2 η|B|u z ûz + a 4 1 2 A 12 • (g z ⊗ ĝ2 z + ĝz ⊗ g 2 z ) -k 2 B 10 (u z ĝz + ûz g z ) + a 5 1 6 A 13 • (g z ⊗ ĝ3 z + ĝz ⊗ g 3 z ) + 1 4 g 2 z : A 22 : ĝ2 z -1 2 k 2 B 20 : (u z ĝ2 z + ûz g 2 z ) -k 2 g z •Q 11 •ĝ z -k 2 ûz 1, X B η B + a 6 1 24 A 14 • (g z ⊗ ĝ4 z + ĝz ⊗ g 4 z ) + 1 12 A 23 • (g 2 z ⊗ ĝ3 z + ĝ2 z ⊗ g 3 z ) + E z •A 11 •ĝ z -1 6 k 2 B 30 • (u z ĝ3 z + ûz g 3 z ) -1 2 k 2 Q 12 • (g z ⊗ ĝ2 z + ĝz ⊗ g 2 z ) -k 2 ûz 1, Y B η B -k 2 u z 1, Y B η B + k 2 u z η|B| D 3 z + O(a 7 ) (52f) (b) Second term of
(52g) J (u; v a ) = a 6 J (u; W ) + o(a 6 )
(c) Third term of (6). The remainder R(u; v a ) in ( 6) can be put in the form

(52h) R(u; v a ) = a,b=R,I ∂Ω 1 0 (1 -t) ∂ ab ϕ u R + tv a R , u I + tv a R , • -∂ ab ϕ u R , u I , • dt v a a v a b dS
The C 0,γ assumption on the derivatives ∂ ab ϕ then yields

(52i) |R(u; v a )| ≤ C v a 2+γ L 2+γ (∂Ω) ≤ C v a 2+γ H 1 (Ω) ≤ Ca 6+3γ = o(a 6
) for some C > 0 and for all a small enough, where we have used the fact that the trace mapping defines a continuous H 1 (Ω) → L p (∂Ω) operator for 1 ≤ p ≤ 4 (e.g. [16, Thm. 6.6-5]), and then invoked Lemma 4.

(d) End of the proof. The theorem follows by using (52f), (52g) and (52i) in ( 6) and grouping all contributions of like orders. Similar arguments would apply, using (50), for components of J defined as integrals over subsets of Ω, if any. 4.2. Centrally symmetric obstacle. Many simple shapes (e.g. balls, ellipsoids or cuboids) fall into this category, for which the following lemma allows expansions to be significantly simpler than (51a-d):

Lemma 5. Let B be centrally symmetric, i.e. x ∈ B =⇒ -x ∈ B, and define the symmetry operator S by Sw(x) = w(-x) for any function w. Then: (i) the following properties hold for any u, w ∈ H 1 (B; C):

(53) (a) SH[u] = H[Su], (d) Su, Sw β B = u, w β B , (b) SN [u] = N [Su],
(e) Su, Sw

η B = u, w η B , (c) SN [u] = N [Su] (ii) if equation (I -H)u B = u holds for some u ∈ H 1 (B; C), then (I -H)Su B = Su also holds.
Proof. To prove (53a), let x ∈ B, implying that -x ∈ B by symmetry assumption. Evaluating SH[u](x) by setting y = -y for the integration variable, we have

SH[u](x) = B β∇G ∞ 0 (-y + x)•∇u B (-y) dV y
Then, property (53a) follows from using in the above equality that (i) [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF]) and (ii) the definition of S implies ∇(Su B )(y) = -∇u B (-y). Similar arguments yield properties (53b,c). Then, properties (53d,e) follow from using the above change of integration variable in definitions [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of misfit function[END_REF]. Finally, item (ii) is obtained at once by applying S to the equation (I -H)u B = u and invoking property (i-a).

∇G ∞ 0 (-y + x) = -∇G ∞ k (y -x) (see
If B is centrally symmetric, we have Sπ m [E m ] = (-1) m π m [E m ] (because π m [E m ] is a homogeneous polynomial of degree m) and SU m B [E m ] = (-1) m U m B [E m ]
(by Lemma 5) for any integer m; in particular, U 1 , U 3 are skew-symmetric while U 2 , U 4 are symmetric. Moreover, S(G1) = G1 by (53b), implying that SX B = X B ; similarly, (G -N )U 1 is skew-symmetric (by (53b,c), and since U 1 is), implying skewsymmetry of Y B -D 3 z . Using these remarks together with properties (53c,d), we find that (54)

A pq = 0, B pq = 0, Q pq = 0 (p + q odd); 1, Y B η B = η|B|D 3 z , 1, Y B η B = η|B| D 3 z
Consequently, T 3 and T 5 are still given by (51a) and (51c), respectively, but (51b) and (51d) reduce to [START_REF] Silva | Higher order topological derivatives in elasticity[END_REF] for the cases p = 1, 2 of relevance here. As a result, the following closed-form expressions for the elastic moment tensors A pq are derived from their definition (33) and using ( 31):

(55) T 4 (z) = 0, T 6 (z) = Re 1 2 J (u; W ) + E z •A 11 •ĝ z -k 2 |B|ηD
(56)

A 11 = β|B| I + βS 1 -1 , A 13 = 3|B| -1 A 11 ⊗ 1, ȳ ⊗ ȳ B , A 22 = β|B| I + βS 2 -1 .
They are then used in formulass (51a), (51c) and ( 55) for the topological derivatives.

4.4. Spherical obstacle. If B is the unit sphere, the volume potentials appearing in the relevant FSTPs (44a), (45a) and in the definition of tensors Q pq can be evaluated in closed form for polynomial densities, see Appendix C; moreover, tensors S 1 , S 2 are given explicitly by (32) while 1, ȳ ⊗ ȳ B = (4π/15)I. As a result, we find that the relevant FSTP solutions are given in B by

U 1 B [E 1 ](x) = 3 β + 3 E 1 •x, U 2 B [E 2 ](x) = 5 2β + 5 E 2 : (x ⊗ x) + β β + 1 1 3 - x 2 2β + 5 E 2 : I, X B (x) = η(2β + 3 -x 2 ) 6(β + 1) k 2 u z
and the associated tensors by

A 11 = 4π β β + 3 I, B 20 = 4π 45 η(5β + 3) β + 1 I, A 13 = 12π 5
β β + 3 I ⊗ I, Q 11 = 4π 5 3η + 2β 2 (β + 3) 2 I, A 22 = 16π 15
β 2β + 5 5I - β β + 1 I ⊗ I , 1, X B η B = 4π 45 η 2 (5β + 6) β + 1 k 2 u z .
5. Discussion.

Computational considerations.

The practical evaluation of the topological derivatives T 3 to T 6 rests on that of the fields u, û, G(•, z), all defined on the same (background) configuration, and of polarization tensors. Moreover, the following remarks apply:

(a) Background and adjoint solutions u, û. They need to be computed only once, irrespective of the number of obstacle sites z considered. They have the representations

(57) u(x) = iρ 0 ω G(• -x) , V D ∂Ω , û(x) = -G(• -x) , ∂ R ϕ -i∂ I ϕ ∂Ω ,
which in practice are explicit only for domains Ω simple enough for G to be known analytically (see Sec. 5.3). As derivatives of u and û of order up to three are involved in T 5 (and orders up to four in T 6 if (54) does not apply), suitable solution or postprocessing methods are needed. If G is known, formulas (57) are well suited for this since they may be differentiated at any order. (b) Polarization tensors. Each tensor needs to be computed only once for given obstacle shape and material properties. This task requires the FSTP solutions

U 1 B , U 2 B .
The latter are known explicitly for ellipsoidal obstacles. For other shapes, the FSTPs (30) for p = 1, 2 must be solved numerically.

(c) Derivatives of the Green's function. Derivatives of either G(•, z) or its complementary part G c (•, z) are needed for the evaluation of T 6 . This warrants closer examination. In most situations, contributions of G(•, z) or G c (•, z) to T 6 have to be computed numerically. This requirement has implications. Firstly, the derivatives ∇ 1 G(z, •) are involved, through (50), in J (u; W ); this looks inconvenient if G is not known analytically as the second argument spans all of ∂Ω. However, the well-known symmetry property

G(z, x) = G(x, z) implies ∇ 1 G(z, x) = ∇ 2 G(x, z),
allowing to evaluate J (u; W ) for given z by means of G(•, z). Using the additive decomposition ( 23) in ( 14), G c (•, z) solves the problem

(a) ∆G c (•, z) + k 2 G c (•, z) = 0 in Ω, (b) ∂ n G c (•, z) = -∂ n G ∞ k (•-z) on ∂Ω.
In fact, upon differentiating that problem w.r.t. the coordinates of z (which is valid since z acts therein as a parameter), the derivatives

H (•, z) := ∇ 2 G c (•, z)•e are found to verify (58) ∆H (•, z)) + k 2 H (•, z)) = 0 in Ω, ∂ n H (•, z) = ∂ n ∂ G ∞ k (• -z) on ∂Ω , 1 ≤ ≤ 3.
which are usual boundary-value problems with smooth boundary data (since z ∈ Ω). Then, T 6 (z) also involves ∇ 12 G c (z, z), which can be evaluated using first-order derivatives of y → H (y, z) at y = z. Secondly, evaluating the expansion of J(a; z) for a large number N of obstacle sites z requires solving many instances of problem (58), which is impractical. A possible remedy consists in using a (truncated

) separable representation of G ∞ k , to have G ∞ k (y-z) = p q=1 α q (y)β q (z) + ε p (
given by e.g. a multipole expansion [START_REF] Gumerov | Fast multipole methods for the Helmholtz equation in three dimensions[END_REF]). The computational work entailed by problems (58) becomes O(p) irrespective of N .

Other boundary conditions on ∂Ω.

Neumann conditions on ∂Ω were assumed for definiteness in Sections 2 to 4, but straightforward adaptations allow to consider other types of boundary conditions. For example, assuming that ∂Ω is divided into complementary subsets S N and S D supporting prescribed values V D of the normal velocity and u D of the acoustic pressure, respectively, the Neumann condition on ∂Ω in problems ( 2) and ( 3) is replaced by ∂ n u = iρ 0 ωV D on S N , u = u D on S D , while the boundary condition (14b) in the definition of G is replaced by

G(•, x) = 0 on S D , ∂ n G(•, x) = 0 on S N .
Due to the modification of the boundary condition setting, the nature of experimental data and the objective function format (3) can also incur changes. For instance, we may here consider the format

J(w) = SN ϕ N w R (y), w I (y), y dS y + SD ϕ D ∂ n w R (y), ∂ n w I (y), y dS y , in (3)
, where ϕ N and ϕ D are C 2 functions with respect to their first two arguments, the new density ϕ D allowing for instance to account for available values V obs of the normal wall acceleration on a measurement surface S obs D ⊂ S D . The derivatives J and J of J are then given by

J (u; w) = p=R,I Sα α=D,N ∂ p ϕ α w p dS, J (u; w) = p,q=R,I α=D,N Sα ∂ pq ϕ α w p w q dS.
instead of (7), while the adjoint field û now satisfies [START_REF] Bendali | Approximation by multipoles of the multiple acoustic scattering by small obstacles in three dimensions and application to the Foldy theory of isotropic scattering[END_REF]. With these changes, Lemma 1 still holds and all results of Secs. 3,4, in particular Theorem 1, are unchanged.

∂ n û = -∂ R ϕ N -i∂ I ϕ N (u R , u I , •) on S N and û = ∂ R ϕ D -i∂ I ϕ D (∂ n u R , ∂ n u I , •) on S D instead of the Neumann condition of problem
5.3. Unbounded media. If an infinite medium is considered (Ω = R 3 ), then of course G = G ∞
k and G c = 0, while the background field u is any field solving (∆ + k 2 )u in R 3 (e.g. a plane wave). This significantly simplifies T 6 since the constants D 3

z and E z appearing in (51d) and ( 55) are then given by D 3 z = η|B|ik 3 u z /4π and E z = -ik 3 g z •A 11 /12π. In that case, all results of Sections 3 and 4 apply with no changes for any objective functional having the format J(w) = D ϕ(w R , w I , •) dV , where D ⊂ R 3 is compact (e.g. D is a measurement region), the adjoint solution û being the acoustic field created by the adjoint source distribution -χ D ∂ R ϕ-i∂ I ϕ , i.e. being given by û

(x) = -G ∞ k (• -x) , ∂ R ϕ -i∂ I ϕ D .
Alternatively, for the halfspace Ω = {y | y 3 < 0} bounded by S = {y | y 3 = 0}, it is well-known that

(59) G c (y, x) = ±G ∞ k (y -x) with x = (x 1 , x 2 , -x 3 )
where the ± sign depends on whether homogeneous Neumann or Dirichlet data is considered on S.

Multiple trial obstacles.

Having so far assumed a single a-dependent obstacle B a , we now briefly explain how Theorem 1 extends to the case of several small obstacles with fixed locations, focusing on the case of two such obstacles (the generalization to three or more being then straightforward). Accordingly, let a second obstacle have support B a := z + aB and relative material parameters β , η ; B a and B a are then scaled by the same characteristic radius a. For this setting, the VIE ( 15) becomes

(60) I -L a u a (x) -L a u a (x) = u(x), x ∈ B a ,
where u a is the restriction of u a to B a and the integral operator

L a : H 1 (B a ) → H 1 (B a ) is defined by L a w (x) = B a η G(•, x)w dV - B a β ∇ 1 G(•, x)•∇w dV
The VIE formulation for the double-obstacle case then consists of equation ( 60) together with its counterpart obtained by reversing the roles of B a and B a . Using an ansatz of the form (35) for each obstacle, the coupling term L a u a (x) has a O(a 3 ) leading contribution, resulting from an outer expansion essentially identical to that of Sec. 3.6. The lowest-order coefficients U 0 , U 1 , U 2 are therefore still given by (39a,b) and (44a) (with similar results for U 1 , U 2 on B ), and are not influenced by the added obstacle B a . The expansion of L a u is in fact formally identical to that of L c a u a , replacing G c (y, x) with G(y, x) and with all other quantities therein now referring to B a , so that we have

P a L a u a (x) = a 3 D 3 z ,z + a 4 D 4 z ,z + E z ,z • x + o(a 4 ) where D 3 z ,z , D 4
z ,z are scalar constants, whose expressions are similar to those of D 3 z , D 4 z , and with

E z ,z = k 2 u z η |B| ∇ 2 G(z , z) -g z •A 11 •∇ 12 G(z , z)
The effect of the coupling on U a thus manifests itself through replacing D 3 z with

D 3 z + D 3 z ,z in (43b) and D 4 z , E z with D 4 z + D 4 z ,z , E z + E z ,z in (43c
) (and the symmetric replacements to obtain governing equations for U 3 , U 4 ). The expansion of objective functionals is insensitive to the new constant D 4 z ,z appearing in U a . By suitable modifications to the derivation of (52c) and (52d), the topological derivatives T p (z, z ) for the dual-obstacle configuration are obtained in terms of the topological derivatives T p (z) and T p (z ) for a single obstacle located at z or z , as given by Theorem 1, as

T 3 (z, z ) = T 3 (z) + T 3 (z ), T 4 (z, z ) = T 4 (z) + T 4 (z ) = 0, T 5 (z, z ) = T 5 (z) + T 5 (z ), T 6 (z, z ) = T 6 (z) + T 6 (z ) + J (u; W, W ) + E z,z •A 11 •ĝ z + E z ,z •A 11 •ĝ z -k 2 |B|ηD 3 z,z û z -k 2 |B |η D 3
z ,z ûz (for centrally-symmetric B, B , and with self-explanatory notation in the expression of T 6 (z, z )). In particular, coupling between obstacles does not occur at orders below O(a 6 ) in the expansion of J(a).

6.

A simple approximate global search procedure. Let J 6 (a; z) be the O(a 6 ) approximation of J(a; z) provided by its high-order topological expansion (Theorem 1), given by (61) J 6 (a; z) := J(0) + a 3 T 3 (z) + a 4 T 4 (z) + a 5 T 5 (z) + a 6 T 6 (z).

We seek the unknown obstacle (assumed to be unique, see discussion in Sec. 5.4) by minimizing the polynomial approximation J 6 (a; z) for sampling points z chosen a priori. This task is, for each z, simple and computationally very light. It can therefore be performed for locations z 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 yielded by this procedure is defined by the location z = x est and size a = R est that minimize J 6 (a; z) for a ≥ 0 and z ∈ G, i.e. given by (62)

x est = arg min z∈G Ĵ6 (z), R est = R(x est ),
with functions z → Ĵ6 (z) and z → R(z) defined through partial minimization of J 6 (a; z) w.r.t. a, i.e.:

(63) Ĵ6 (z) = min a≥0 J 6 (a; z), R(z) = arg min a≥0 J 6 (a; z).

In all numerical examples to follow, the "asymptotic surrogate" J 6 of J is set up on the assumption of a spherical trial scatterer B a . We note in addition that a similar procedure is developed and implemented for 2D EIT problems in [START_REF] Ferreira | A new non-iterative reconstruction method for the electrical impedance tomography problem[END_REF]. The true obstacle B (Fig. 1) is centered at x = (2.05d, 1.25d, -2.05d), its material contrasts being β = 1, η = -0.5. Three geometries are considered for B: a sphere (S) of radius d/5, a horizontally elongated ellipsoid (E) whose semiaxes (aligned with the (x 1 , x 2 , x 3 ) axes) are (A, 1 2 A, 1 2 A) with A = 2 2/3 d/5, and a banana shape (B) obtained by applying the transformation x 3 -→ x 3 -0.15d × (x 1 -x1 ) 2 /A 2 to (E). For ease of comparison, the 'true' radius R is defined as that of the sphere having the same volume as B, so that R = d/5 for the three cases. Three frequencies are considered, such that kd = 1, 2, 5. Note that x ∈ G: four sampling points z ∈ G are closest to x, with |z -x| = √ 5d/20 ≈ .1118d. Except in the last paragraph (e) of this section, the relative contrasts of the trial obstacle B a are β = β, η = η.

(a) Obstacle size estimation (known location, noise-free data). In this preliminary numerical experiment, the size of an obstacle whose location is known is estimated by computing R(x) defined by (63). Results for the relative error R(x)/ R -1 on the obstacle radius found using noise-free synthetic data are given in Table 1 for all of the above-defined scatterer configurations, testing configurations and frequencies. The size estimation accuracy is seen to decrease as the frequency increases, and to be relatively insensitive to the density of the measurement and testing grids.

(b) Approximate global search, noise-free data. The approximate global search procedure outlined in Section 6 has then been performed on a search grid G of N G = 51 × 51 × 51 regularly spaced sampling points spanning the 3-D box-shaped search region S defined by -5d ≤ x 1 , x 2 ≤ 5d, -5.5d ≤ x 3 ≤ -0.5d. The polynomial approximation J 6 (a; z) of J has been set up for all N G sampling points of G and using the explicit Green's function.

The estimated obstacle radius R est defined by (62) obtained for all obstacle geometries, testing configurations and frequencies and using noise-free synthetic data is compared to R in Table 2. The size estimation accuracy again decreases as the frequency increases, and is relatively insensitive to the density of sources and sensors: even the 2 × 2 testing configuration, which has only 4 sources and 9 sensors, yields good results. For all cases displayed in Table 2, the identified obstacle location x est is closest to x among the available grid points. For cases (E) and (B) involving 'true' scatterer shapes that increasingly deviate from the trial spherical shape, the accuracy for the 'equivalent radius' R (i.e for the obstacle volume) is nonetheless similar to that achieved for case (S).

Additionally, values of Ĵ6 (z) close to the minimum J min 6 := Ĵ6 (x est ) are found to occur only at grid points close to x est (with R(z) ≈ R est at those locations). This is illustrated in Figure 2, where iso-surfaces Ĵ6 (z) = ζJ min 6 , with ζ = 0.6 , 0.7 , 0.8 , 0.9, depicted for obstacle configuration (E) and testing configuration 5 × 5, are seen to shrink to a small neighbourhood of x as Ĵ6 approaches J min The measurement residuals u B e -u obs e , which carry the useful information for obstacle identification, are severely affected by even small perturbations of u obs e since their magnitude is much smaller than that of the measurement u obs , especially at the lower frequency kd = 1. Estimations R est and x est have been computed for cases (S), (E) and (B), using the same sampling grid G, for all testing configurations and frequencies. The relative errors on R est and offsets |x est -x| achieved are shown in Table 3 for noise level σ = 0.02. On comparing these results with those for errorfree data (Table 2), the accuracy deterioration remains moderate in most cases with kd = 2 or kd = 5, while being unacceptable for the lower frequency kd = 1. In the former case, the scarcest testing configuration 2 × 2 also no longer yields satisfactory results. These trends are, unsurprisingly, more pronounced still for stronger data corruption σ = 0.05 (Table 4), where reasonable values are obtained only for kd = 2 and kd = 5 with testing configurations 5×5 or 10×10. On the other hand, identification results given in Table 5 for case (S) using synthetic data with 20% relative noise on scattered field are quite similar to those obtained using noisefree data (Table 3), suggesting robustness of the search method against significant perturbation of the essential data. Table 3. Offset |x est -x| and relative error ε R := R est / R -1 on obstacle radius for obstacles (S), (E) and (B): testing configurations 2 × 2, 5 × 5 and 10 × 10, synthetic data with 2% relative noise on total field. Where |x est -x| is unacceptably large, ε R is deemed irrelevant and not shown. (d) High-order expansion vs. topological derivative. As seen in previous examples, the accuracy of the obstacle size estimation using the high-order expansion J 6 decreases as the probing frequency is increased. On the other hand, the qualitative imaging provided by the usual topological derivative T 3 is known to deteriorate as the probing frequency decreases, see e.g. [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]. These differing trends are further exemplified in Figures 3 and4, which show contour plots of mappings z → T 3 (z) and z → Ĵ6 (z) treated as imaging functions. The image given by Ĵ6 (z) at low frequencies is less smeared and better pinpoints the true location of the unknown object. Both images get sharper as frequency increases.

ka = 1 ka = 2 ka = 5 |xest -x| ε R |xest -x| ε R |xest -x| ε R 2 × 2 6.5e+00 -- 4
(e) Varying the physical parameters of the trial obstacle. The higher-order topological expansion J 6 depends on the assumed shape B and material parameters ρ , c of the trial obstacle B a . This dependence has not been exploited so far, and in-depth analysis of its effect and potential benefits on the proposed search procedure is left to future work. As a preliminary illustration of whether that procedure is sensitive to the choice of trial material parameters, Figure 5 shows a contour plot of (β, η) → Ĵ6 (x; β, η) for obstacle configuration (B), testing configuration 5 × 5, probing frequency kd = 2 and noise-free data (the true material parameters being ( β, η) = (1, -0.5)). A significant proportion of the possible choices of (β, η) (the dark blue region of the contour plot) were found to be such that Ĵ6 (x; β, η) = J 6 (0; x) ≈ 3.6 10 -6 (i.e. a → J 6 (a; x) has no strictly positive minimizer). The sensitivity of Ĵ6 (x; β, η) to (β, η) in the vicinity of ( β, η) is moreover seen to be rather mild, which suggests that accurate simultaneous determination of the obstacle size and material parameters might not be easy.

8. Proof of Proposition 1, part 2. Regarding item (i), the boundedness and invertibility of A a result from part 1 and the definition (18) of the scaling operator P a . To estimate A -1 a , we consider the equation A a u a = u for some given background field u ∈ H 1 (B a ). This equation is equivalent to (I -H)U a = P a [u] with U a := P a [u a ]. It can be solved for U a in two steps. Step (i) consists in solving for ∇U a the singular VIE ∇U a -∇HU a = ∇P a [u], which is known [START_REF] Bonnet | A modified volume integral equation for anisotropic elastic or conducting inhomogeneities. Unconditional solvability by Neumann series[END_REF] to be for any a ≤ a 0 .

Summing inequalities (65a) and (65b), we have that ∇u a H 1 (Ba) ≤ C 0 ∇u H 1 (Ba) for some C 0 > 0, which completes the proof of item (i).

Regarding item (ii), the integral operator L c a (given by ( 25)) has kernels that are either bounded or weakly singular over B × B; it is therefore a bounded (in fact, compact) H 1 (B) → H 1 (B) operator by e.g. [START_REF] Hsiao | Boundary Integral Equations[END_REF]Thm. 6.1.12]. Moreover, these kernels being O(1) as a → 0 (see the kernel expansions given after (39) in Sec. 3.5.1), there exists C c > 0 such that L c a ≤ C c for a small enough, as claimed. Then, using equalities (47), we have (66)

L c a u H 1 (Ba) = a 2 P -1 a L c a P a u H 1 (Ba) ≤ a 5/2 L c a P a u H 1 (B) ≤ a 5/2 C c P a u H 1 (B) ≤ aC c u H 1 (Ba)

Therefore, there exists C c such that L c a ≤ C c for any a < a 0 . The proof of item (ii) is complete. 9. Proof of Proposition 2. As a consequence of equation ( 15) satisfied by u a , the inner expansion error δ a := u a -[P -1 a U Proof of Lemma 6. The decomposition I -L a = A a -L c a holds, with A a and L c a as in Proposition 1. Since A a is invertible (part 2(i) of Prop. 1), we write I -L a = A a I-A -1 a L c a . Now, since L c a ≤ aC c for any a < a 0 (part 2(ii) of Prop. 1), there exists for any c > 1 an inhomogeneity size a 1 such that A -1 a L c a H 1 (Ba) ≤ c < 1 for any a ≤ a 1 , namely a 1 = min a 0 , (C 0 C c ) -1 c with C 0 as in Prop 1. For all a < a 1 , I -A -1 a L c a is then invertible by Neumann series, so has a bounded inverse:

I -A -1 a L c a -1 H 1 (Ba) ≤ 1 -A -1 a L c a H 1 (Ba) -1 ≤ (1 -c) -1
Concluding, I -L a is invertible and its inverse, given by (I -L a ) -1 = I -

A -1 a L c a -1 A -1
a , is bounded uniformly in a for all a < a 1 .

Proof of Lemma 7. We first derive an explicit expression of γ a , given by (67), recalling the definition (35) of U

a . First, using the rescaled form [START_REF] Laurain | Topological sensitivity analysis in fluorescence optical tomography[END_REF] of L a , we have P a (I -L a )[P -1 a U (4) a ] = I -H)U (4) a -a 2 L c a U (4) a . Then, effecting the combination (38a)+a(38b)+a 2 (43a)+a 3 (43b)+a 4 (43c) of the VIEs governing the coefficients U 0 , . . . , U 4 , I -H)U (4) a is found to be given by (I -H)U (4) a = T 4 [u] + a 2 L c 0 U 0 + a 3 (L c 0 U 1 + L c 1 U 0 ) + a 4 (L c 0 U 2 + L c 1 U 1 + L c 2 U 0 ) with the Taylor polynomial T 4 [u] as defined in (37). Substituting this into the above expression of P a (I -L a )[P -1 a U

a ], using the resulting formula for evaluating P a γ a with γ a as given by (67), and rearranging terms, we obtain Using estimates (69a), (69b) and (69c) in (68) therefore yields P a γ a H 1 (B) ≤ a 5 C for some C > 0. Lemma 7 finally follows from applying relations (47) to the latter estimate of P a γ a .

  )). Moreover, since L c a = O(1) (see Prop. 1), it is convenient to proceed as follows: (i) state and solve the O(1) and O(a) equations, which are simple and yield U 0 and U 1 ; (ii) compute the expansion of L c a U (4) a and use therein the known values of U 0 , U 1 ; (iii) set up and solve for U 2 , U 3 , U 4 the more-involved O(a 2 ) to O(a 4 ) equations. (i) Determination of U 0 , U 1 . The O(1) and O(a) equations arising from (36) are (38) (a) I -H U 0 = u z , (b)

Lemma 4 .

 4 Let D Ω be an open subset of Ω, independent of a, such that z ∈ D, and assume a is small enough to have B a D. Then: v a -a 3 W H 1 (Ω\D) = O(a 4 ). Proof. Let D D be an open subset of D such that B a D , and define w := (v a -a 3 W )χ in Ω (with χ a smooth cut-off function, equal to 1 in Ω \ D and vanishing in D ). Let b := -(∆ + k 2 )w. Since (∆ + k 2 )v a = (∆ + k 2 )W = 0 in Ω -D, supp(b) ⊂ D \ D . Expansion (50) hold pointwise in D \ D , and so does the corresponding expansion of ∇v a by virtue of y → ∇ 2 G(y, x) being smooth in B a , implying that b L 2 (Ω) = b L 2 (D\D ) = O(a 4 ). As w solves -(∆ + k 2 )w = b in Ω and ∂ n w = 0 on ∂Ω, it depends linearly and continuously on b, so that w H 1 (Ω\D) = O(a 4 ).

  straightforwardly. In (52b) and for the remainder of this proof, we use for convenience the shorthand notationsU m B := U m B [g m z ]

7 .

 7 Numerical experiments. To demonstrate the above approximate search procedure, we consider the identification of a penetrable object ( B, ρ,c) embedded in an acoustic medium occupying the half-space Ω = {x |x 3 ≤ 0} (Fig.1). A homogeneous Neumann condition is assumed on the surface S = {x | x 3 = 0}. The relevant Green's function G is then given in closed form by[START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF] and (59).Three synthetic limited-aperture testing configurations (labelled 2 × 2, 5 × 5 and 10 × 10 in the sequel) are defined, where the squareregion {x | -5d ≤ x 1 , x 2 ≤ 5d, x 3 = 0} of S is divided into 2 × 2,5 × 5 and 10 × 10, respectively (d being an arbitrary reference length). 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 and 10 × 10 feature N = 4, 25, 100 sources and M = 9, 36, 121 sensors, respectively. A set of N synthetic experiments is defined, each consisting in activating one of the N sources. The identification is formulated in terms of the cumulative least-squares cost functional (64) J (B, β, η) m ) -u obs e (x m ) 2 where u obs e and u B e are the acoustic fields induced by point source x e for the 'true' and 'trial' objects ( B, β, η) and (B, β, η). The polynomial approximation J 6 (a) of J(a) := J (B a , β, η) defined by (61) has been set up for the functional (64) and a spherical trial scatterer, i.e. using (55) and the formulas of Sec. 4.4 to evaluate the coefficients T 3 (x), . . . , T 6 (x) given by Theorem 1. The background and adjoint fields for the e-th experiment are given by u e = G(•, x e ), ûe = M m=1 u B e (x m ) -u obs e (x m ) G(•, x m )where the overbar indicates complex conjugation. The synthetic data u obs e are computed by means of a coupled system of boundary integral equations (BIEs), the boundary of B being meshed using 600 eight-noded boundary elements, leading to 3604 nodal unknowns for the BIE system.

Figure 1 .

 1 Figure 1. Identification of a penetrable scatterer ( B, ρ,c) in a acoustic half-space: geometry and notation (the dark shaded part is the search region S).

Table 2 . 5 - 2 .

 252 Relative error R(x)/ R -1 on obstacle radius for obstacles (S), (E) and (B) of unknown location: testing configurations 5 × 5, 10 × 10 and 20 × 20, noise-free synthetic data. A distance |x est -x| = d √ 5/20 is found for all cases. 1e-02 -6.0e-02 -2.4e-01 10 × 10 -2.2e-02 -6.2e-02 -2.5e-01 (E) 2 × 2 -4.7e-03 -4.6e-02 -2.3e-01 5 × 5 -3.5e-03 -5.2e-02 -2.4e-01 10 × 10 -3.4e-03 -3.5e-02 -2.4e-01 (B) 2 × 2 -2.2e-03 -3.5e-02 -2.1e-01 5 × 5 -1.2e-03 -3.7e-02 -2.2e-01 10 × 10 -1.1e-03 -3.9e-02 -2.2e-01

Figure 2 . 6 ,

 26 Figure 2. Iso-surfaces of Ĵ6 (z) for Ĵ6 = ζJ min 6 , with ζ = 0.6 (top left), ζ = 0.7 (top right), ζ = 0.8 (bottom left) and ζ = 0.9 (bottom right): obstacle (E), testing configuration 5 × 5 and noise-free data. The iso-surfaces and true obstacle location are emphasized by projections on coordinate planes.

6 .

 6 (c) Approximate global search, noisy data. Then, the influence of data noise on the search procedure is considered, by replacing the (synthetically) measured total field u obs e entering the cost functional (64) 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 ) , η e,m , η e,m being uniform random numbers with zero mean and standard deviation σ.

Figure 3 .

 3 Figure 3. Contour plots of z → T 3 (z) (left column) and z → Ĵ6 (z) (right column) in the horizontal plane containing the true obstacle center x, for obstacle (B), testing configuration 5 × 5 and noise-free data. Testing frequencies are kd = 1 (top row), kd = 2 (middle row) and kd = 5 (bottom row).

Figure 4 .

 4 Figure 4. Contour plots of z → T 3 (z) (left column) and z → Ĵ6 (z) (right column) in a vertical plane containing the true obstacle center x, for obstacle (B), testing configuration 5×5 and noise-free data. Testing frequencies are kd = 1 (top row), kd = 2 (middle row) and kd = 5 (bottom row).

Figure 5 .

 5 Figure 5. Sensitivity of search procedure to trial physical parameters: contour map of (β, η) → Ĵ6 (x; β, η) (actual obstacle parameters are β = 1, η = -0.5)

Lemma 6 .Lemma 7 .

 67 )a ] ∈ H 1 (B a ) solves the integral equation (67) (I -L a )δ a = γ a , with γ a := u -(I -L a )[P -1 a U a ] Proposition 2 follows directly from the following two lemmas. Let β such that -1 < β < ∞. The integro-differential operator I -L a : H 1 (B a ) → H 1 (B a ) is invertible. Moreover, there exists a 1 > 0 such that (I -L a ) -1is bounded uniformly in a for all a ≤ a 1 . The right-hand side γ a of equation (67) satisfies γ a H 1 (Ba) = O(a 11/2 ).

P

  a γ a = P a u -T 4 [u] (68) + a 2 L c a -L c 0 -aL c 1 -a 2 L c 2 U 0 + a 3 L c a -L c 0 -aL c 1 U 1 + a 4 (L c a -L c 0 )U 2 + a 5 L c a (U 3 + aU 4 ) We now estimate P a γ a L 2 (B)and ∇(P a γ a ) L 2 (B) . Firstly, the kernel functions of the integral operatorsL c a -L c 0 , L c a -L c 0 -aL c 1 and L c a -L c 0 -aL c 1 -a 2 L c2 can be shown to be bounded for (ȳ, x) ∈ B × B and (using Taylor expansions in a about a = 0) to have respective leading orders O(a), O(a 2 ) and O(a 3 ). Therefore, they are continuous H 1 (B) → H 1 (B) operators with O(a 2 ) continuity constants for a < a 1 , and we havea 2 L c a -L c 0 -aL c 1 -a 2 L c 2 U 0 + a 3 L c a -L c 0 -aL c 1 U 1 (69a) + a 4 (L c a -L c 0 )U 2 H 1 (B)≤ a 5 C for some C > 0. Then, since U 3 , U 4 solve FSTPs whose definition does not depend on a and whose data has H 1 (B) regularity while L c a is (by Prop. 1) a bounded H 1 (B) → H 1 (B) operator there exists a constant C > 0 depending only on B and β such that, for any a < a 1(69b) L c a (U 3 + aU 4 ) H 1 (B) ≤ C.Finally, provided u has sufficient (at least C 5 ) local regularity at z, we have (69c) P a u -T 4 [u] H 1 (B) = O(a 5 )

Table 1 .

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

		kd = 1	kd = 2	kd = 5
	(S) 2 × 2	-1.1e-02 -4.7e-02 -2.1e-01
	5 × 5	-1.1e-02 -4.9e-02 -2.2e-01
	10 × 10 -1.1e-02 -5.0e-02 -2.3e-01
	(E) 2 × 2	-1.8e-03 -3.7e-02 -2.1e-01
	5 × 5	-5.0e-03 -4.1e-02 -2.2e-01
	10 × 10 -5.8e-03 -4.2e-02 -2.2e-01
	(B) 2 × 2	-4.4e-03 -4.0e-02 -2.1e-01
	5 × 5	-8.2e-03 -4.4e-02 -2.2e-01
	10 × 10 -9.2e-03 -4.5e-02 -2.3e-01

Table 4 .

 4 Offset |x est -x| and relative error ε R := R est / R -1 on obstacle radius for obstacles (S), (E) and (B): testing configurations 2 × 2, 5 × 5 and 10 × 10, synthetic data with 5% relative noise on total field. Where ε R is unacceptably large, R est / R -1 is deemed irrelevant and not shown.

			ka = 1		ka = 2		ka = 5
			|xest -x|	ε R	|xest -x|	ε R	|xest -x|	ε R
		2 × 2	3.5e+00	--	6.9e+00	--	3.9e+00 -2.4e-01
	(S)	5 × 5	3.7e+00	--	7.2e+00	--	1.1e-01 -2.3e-01
		10 × 10	5.6e+00	--	8.5e+00	--	1.1e-01
		2 × 2	5.1e+00	--	5.0e+00	--	4.2e+00	--
	(E)	5 × 5	6.4e+00	--	3.7e+00	--	1.1e-01 -1.9e-01
		10 × 10	4.9e+00	--	5.4e+00	--	1.1e-01 -2.1e-01
		2 × 2	5.9e+00	--	3.8e+00	--	6.1e+00	--
	(B)	5 × 5	4.8e+00	--	5.1e+00	--	1.1e-01 -2.1e-01
		10 × 10	3.8e+00	--	6.0e+00	--	1.1e-01 -2.1e-01

Table 5 .

 5 Offset |x est -x| and relative error ε

		ka = 1		ka = 2		ka = 5
		|xest -x|	ε R	|xest -x|	ε R	|xest -x|	ε R
	2 × 2	1.1e-01	6.7e-03	1.1e-01 -5.3e-02	1.1e-01 -2.3e-01
	5 × 5	1.1e-01 -2.1e-02	1.1e-01 -6.6e-02	1.1e-01 -2.4e-01
	10 × 10	1.1e-01 -2.1e-02	1.1e-01 -6.3e-02	1.1e-01 -2.5e-01

R := R est / R -1 for obstacle (S) of unknown location: testing configurations 2×2, 5×5 and 10 × 10, synthetic data with 20% relative noise on scattered field.

Appendix A. Concise derivation of VIE [START_REF] Cedio-Fengya | Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction[END_REF]. Testing equation (14a) by any function v ∈ H 1 (Ω) ∩ C 1 (D x ), where D x is a neighbourhood of x, and applying the first Green's identity to the resulting integral over Ω yields the identity (I):

v, G(•, x) Ω -k 2 v, G(•, x) Ω = w(x) (the domain integrals being well defined since ∇G(•, x) ∈ C ∞ (Ω \ {x}; C 3 ) ∩ L 1 (Ω; C 3 )). Then, noting that both u and u a have the requisite local regularity at x for any x ∈ B a ∪ (Ω \ B a ), we write problems [START_REF] Ammari | Polarization and Moment Tensors with Applications to Inverse Problems and Effective Medium Theory[END_REF] with B = B a and (2) in weak form (both with test function G(•, x)) and use (I) with w = u a and w = u, respectively, to obtain

Then, the VIE ( 15) is obtained by subtracting (b) from (a) and expressing the integrals over B a in terms of the operator L a introduced in (16).

Appendix B. FSTP solutions for ellipsoidal inhomogeneities. Let B be an ellipsoid. Explicit derivations (by methods given in e.g. [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] for elastic inclusions) then show that if U is polynomial with degree m in B and vanishes outside of B, the restriction to B of HU is polynomial with degree m. In particular, for densities

where S 1 (B) and S 2 (B) are constant tensors (respectively of order 2 and 4) for which analytical expressions are available, known for the case of elastic inhomogeneities as Eshelby tensors; moreover, they depend only on the shape (i.e. aspect ratios) of B, not on its size (i.e. S i (λB) = S i (B) for i = 1, 2 and any λ > 0). If B is a ball, S 1 and S 2 are given by (32), as shown in Appendix C. It is then a simple matter to check, using (70), that the solution [START_REF] Silva | Higher order topological derivatives in elasticity[END_REF].

Appendix C. Analytical evaluation of volume potentials. Exact solutions of FSTPs for the unit ball B can be derived from analytical expressions of the volume potential U[g], defined by

for polynomial densities g. First, setting y = x + r and R(z

with r := |r|, r = r/r, and Ŝ being the unit sphere in R 3 . Let g 0 (y) = 1, g 1 (y) = y, g 2 (y) = y ⊗ y. Then:

Elementary integrations and rearrangements yield (with

and the same treatment is applicable to polynomial densities g of higher degree. Then, the volume integral operator H evaluates as

, we obtain Appendix D. Auxiliary proofs. Operators G, H and M, N verify the following interrelations and symmetry properties, which will be repeatedly used:

We use in the sequel the shorthand notations introduced in the proof of Theorem 1, namely

By applying the operator β∇ to the VIE [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Equations[END_REF], taking the inner product of the resulting equality with ∇u and integrating over B, we find u , u B -Hu

The Lemma then follows by writing the same identity with the roles of (u, u B ) and (u , u B ) reversed, subtracting both equalities and invoking property (72c) of H. D.2. Proof of Lemma 3. First identity:

Second identity:

by Lemma 2 and (45b)

Third identity: We have

where the second equality follows from identities 72) and the third uses definitions (46a) of B 21 and (46b) of Q 21 . Now, expressing X B by means of equation (45a) and using the interrelation (72a), we have

so that the last three terms inside the square brackets of (73) are given by Using the above equality in (73) finally yields the desired identity.