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

. This article concerns an extension of the topological derivative concept for 3D inverse acoustic scattering problems involving the identiﬁca- tion of penetrable obstacles, whereby the featured data-misﬁt cost function J is expanded in powers of the characteristic radius a of a single small inhomo- geneity. The O ( a 6 ) approximation J 6 of J is derived and justiﬁed for a single obstacle of given location, shape and material properties embedded in a 3D acoustic medium of arbitrary shape. The generalization of J 6 to multiple small obstacles is outlined. Simpler and more explicit expressions of J 6 are obtained when the scatterer is centrally-symmetric or spherical. An approximate and computationally light global search procedure, where the location and size of the unknown object are estimated by minimizing J 6 over a search grid, is proposed and demonstrated on numerical experiments, where the identiﬁcation from known acoustic pressure on the surface of a penetrable scatterer embedded in a acoustic semi-inﬁnite medium, and whose shape may diﬀer from that of the trial obstacle assumed in the expansion of J , is considered.


1.
Introduction. Inverse scattering has been intensively investigated over the last quarter century, in particular due to the development of qualitative, sampling-based, methods [14,23,29] 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. [5,6,7,13,20,21,24,26]. 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 [17,20], 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 [30]). 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. [2,8], 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 [9] for sound-hard obstacles in 3D acoustic media, [10] for cracks in 2D conducting media, [18] for 2D electrical impedance tomography (EIT). Elastic media are considered in [31] (expansion of the potential energy in 2D solids with small holes) and [12] (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 [9,10] 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 [9,10], 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 (1) 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 (1). 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 [9,10,31] 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 [9] for sound-hard obstacles, and also used in [18] 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 [9,10], 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 [2] using coupled boundary integral equation formulations. Matched expansions are used in [4] (to obtain uniform O(a 3 ) expansions for penetrable objects) and [8] (who give arbitrary-order expansions for impenetrable objects). Arbitrary-order expansions are given in [3] 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 [11,25], 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 [12]. 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.
2. 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 Ω, 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 (2) 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): 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 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 (4); 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) 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 (4). Denoting by u a (·; z) the total field solving problem (3), we define the cost function J(a; z) in terms of J by (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.: Lemma 1. Let u a andû be the total and adjoint fields, respectively solving problems (3) with B = B a and (9). The first-order derivative J (u; v a ) in (7) is then given by Proof. Writing in weak form the adjoint problem (9) (with test function v a ) and the forward problems (2) and (3) (with test functionû), we have the identities 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 [17,20] as 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. [3,20,1] (see e.g. equation 57 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 where U 1 , W are known functions that depend on B and β, η (see equations (40b) and (51)).
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 (7) and (13), the lowest expansion order meeting this goal is O(a 6 ). Moreover, the known outer solution expansion (13) suffices for determining the (leading) O(a 6 ) expansion of J (u; v a ) through (8), 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 (12). 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 (18) implies ∇u a (x) = a −1 ∇U a (x), we obtain (14) 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 (14) is consistent with the expected leading-order expansion (11).
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 (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) in which I is the identity operator and the integral operator L a is defined by 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 (16) is given in Appendix A; otherwise see [11] and references therein. In the case Ω = R 3 , for which G = G ∞ k (see eq. (24)), the VIE (16) is identical (up to notational adjustment) to that established in [25] and applied to the obstacle configuration (B a , β, η).

3.2.
Rescaled form of VIE formulation. The integral equation (16) is now recast in rescaled form (i.e. using integral operators defined on the fixed compact domain B) by effecting the coordinate change (18) x noting that this process also transforms the differential volume element according to dV x = a 3 dVx (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 (18) to functions, and its inverse P −1 a : . The coordinate stretching exerted by (18) then implies that gradients transform under P a according to Applying the scaling (18) to both arguments x, y in the VIE (16), we obtain the governing VIE 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 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 , the singular fundamental solutions of the Helmholtz and Laplace equations, given by The complementary acoustic Green's function G c (·, x) is seen (using superposition) to solve a BVP of the form (2) with boundary condition , 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 (18) implies . The above properties, combined with (24), allow to recast the scaled integral operator L a given by (23) as (25) L a = H + a 2 L c a , where the leading operator H and the complementary operator L c a are defined by noting that H does not depend on a. The VIE (21) for U a then takes the form 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: Let a 0 such that B a Ω for any a ≤ a 0 . Then: Regarding item 1, (a) stems from Thm. 6.1.12 in [22] on pseudodifferential operators, and (b) is [11,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 (18) on density gradients, see (20).
The rescaled version (27) of the VIE (16) 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 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. [11]): given a background field u ∈ H 1 loc (R 3 ), find the total field u B such that Solutions to such FSTPs satisfy the following reciprocity identity (see proof in Appendix D.1): . 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 (30) u 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 (30) 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 The above definitions, and in particular this convention, imply , which satisfies the well-posed integral equation . We will need cases p = 1, 2, for which (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) Polarization tensors (PTs). They appear in many asymptotic expansions involving small inhomogeneities, see e.g. [3,15]. For integers p, q ≥ 1, let the PT A pq be the constant tensor of order p + q such that for any constant tensors E p and E q of respective orders p, q (equation (34) 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 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 [3] (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 (22), the desired inner expansion of u a in B a is assumed to have the form 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 (36) into the asymptotic expansion form (27) of the VIE, i.e. by writing that the equation 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 (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 (30) with g z := g 1 z . Having the form (31), respectively with p = 0 and p = 1, their solutions are (see Sec. 3.4) We note in particular that ∇U 0 = 0, which will bring some simplification in forthcoming expansions and shows that the leading contribution of (14) is, as expected (see Sec. 2.5), of order O(a 3 ).
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 (26) We note that all terms in formulas (41) except those involving operators M, N , G yield constant or linear functions ofx for any U . The above definitions then result in the following expansion of L c a U , wherein the scalar constants D 3 z , D 4 z ∈ C and the constant vector E z = C 3 can be evaluated by using (40) and expressing relevant integrals on B in terms of the polarization tensors defined by (34) and (47a) whenever possible. Introducing for convenience the shorthand notations G c whereas the precise value of D 4 z will prove irrelevant.
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 where N is the operator defined in (42). We note that B 01 = 0 (due to U 0 B being a constant function and z being the centroid of B by assumption).
with the tensors B pq , Q pq as defined by (47a,b) and where D 3 z and Y B (in the last identity) are defined by (43) and problem (46b) with u z , g z replaced withû z ,ĝ z .

3.5.2.
Resulting solution expansion and its justification. To assess the approximation order achieved by (36), 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 (19)) implies 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 (36), 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 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 (16) used as an integral representation: Since x ∈ B a , y → G(y, x) is smooth for y ∈ B a , and one therefore has 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. [2,4]. Expansion (51) 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. [27,Thm. 4.18]). Such extra smoothness for ∂Ω is avoided by the following version of the outer expansion: 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 (6) 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 with the topological derivatives T 3 , . . . T 6 given by The tensors A pq , B pq and Q pq are respectively defined by (34), (47a) and (47b) in terms of solutions to free-space transmission problems (FSTPs) with polynomial background field (see Sec. 3.4), the function W is given by (51), the functions X B , Y B solve the FSTPs (46a,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 (38). Proof. The proof results from expanding in powers of a each term of expansion (7) of J(u a ).
(a) First term of (7). Recalling Lemma 1 and eqn. (14), and effecting the coordinate change (18) in the relevant integrals, J (u; v a ) is given by of U a is given by . We first evaluate the term P aû , U (4) a β B . Using (53b) and a Taylor expansion in a of P aû (i.e. (38) with u replaced withû), introducing the polarization tensors A pq and invoking the symmetry property (35) wherever possible, we obtain 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 : , 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 Similarly evaluating the term P aû , U (4) a η B by using the inner expansion (53b), 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 Moreover, the second term in the right-hand side of (53a) admits the following estimate by virtue of Proposition 2 and the Cauchy-Schwarz inequality: (7). It admits, as a direct consequence of the outer expansion (51), the expansion (c) Third term of (7). The remainder R(u; v a ) in (7) can be put in the form The C 0,γ assumption on the derivatives ∂ ab ϕ then yields ) 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 (53f), (53g) and (53i) in (7) and grouping all contributions of like orders. Similar arguments would apply, using (51), 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 (52a-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): Proof. To prove (54a), let x ∈ B, implying that −x ∈ B by symmetry assumption. Evaluating SH[u](x) by setting y = −y for the integration variable, we have Then, property (54a) follows from using in the above equality that (i) ∇G ∞ 0 (−y + x) = −∇G ∞ k (y − x) (see (24)) and (ii) the definition of S implies ∇(Su B )(y) = −∇u B (−y). Similar arguments yield properties (54b,c). Then, properties (54d,e) follow from using the above change of integration variable in definitions (10).  (by (54b,c), and since U 1 is), implying skewsymmetry of Y B − D 3 z . Using these remarks together with properties (54c,d), we find that (55) Consequently, T 3 and T 5 are still given by (52a) and (52c), respectively, but (52b) and (52d) reduce to Ellipsoidal obstacle. Ellipsoids are centrally-symmetric shapes for which ∇U p B [E p ] are polynomials, see (32) 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 (34) and using (32): They are then used in formulass (52a), (52c) and (56) for the topological derivatives.

Spherical obstacle.
If B is the unit sphere, the volume potentials appearing in the relevant FSTPs (45a), (46a) 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 (33) while 1,ȳ ⊗ȳ B = (4π/15)I. As a result, we find that the relevant FSTP solutions are given in B by

and the associated tensors by
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 , 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 (55) does not apply), suitable solution or postprocessing methods are needed. If G is known, formulas (58) 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 (31) 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 (51), 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, allowing to evaluate J (u; W ) for given z by means of G(·, z). Using the additive decomposition (24) in (15), 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 (59) ∆H (·, z)) + k 2 H (·, z)) = 0 in Ω, 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 (59), 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 [19]). The computational work entailed by problems (59) becomes O(p) irrespective of N .

5.2.
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 (15b) 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 (4) can also incur changes. For instance, we may here consider the format in (4), 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 instead of (8), while the adjoint fieldû now satisfies on S D instead of the Neumann condition of problem (9). With these changes, Lemma 1 still holds and all results of Secs. 3,4, in particular Theorem 1, are unchanged.

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 (52d) and (56) 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û( Alternatively, for the halfspace Ω = {y | y 3 < 0} bounded by S = {y | y 3 = 0}, it is well-known that 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 (16) becomes where u a is the restriction of u a to B a and the integral operator L a : The VIE formulation for the double-obstacle case then consists of equation (61) together with its counterpart obtained by reversing the roles of B a and B a . Using an ansatz of the form (36) 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 (40a,b) and (45a) (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 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 [18]. 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 (24) and (60).
Three synthetic limited-aperture testing configurations (labelled 2 × 2, 5 × 5 and 10 × 10 in the sequel) are defined, where the square region {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  a penetrable scatterer (B,ρ,c) in a acoustic half-space: geometry and notation (the dark shaded part is the search region S). Table 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. boundary ofB being meshed using 600 eight-noded boundary elements, leading to 3604 nodal unknowns for the BIE system.
(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 (64). Results for the relative error R(x)/R − 1 on Table 2. 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.  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 (63) obtained for all obstacle geometries, testing configurations and frequencies and using noise-free synthetic data is compared toR 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 tox 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 ofx asĴ 6 approaches J min 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 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.  Table 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.
3.8e+00 --6.0e+00 --1.1e−01 −2.1e−01 Table 5. Offset |x est −x| and relative error ε 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. (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. [17,20]. These differing trends are further exemplified in Figures 3 and 4, 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.
(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 Figure 3. Contour plots of z → T 3 (z) (left column) and z → J 6 (z) (right column) in the horizontal plane containing the true obstacle centerx, 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). J 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.   Step (ii) then uses the representation U a = HU a +P a [u] in B, yielding U a L 2 (B) ≤ P a u L 2 (B) + C ∇(P a u) L 2 (B) and thus u a L 2 (Ba) ≤ u L 2 (Ba) + Ca ∇u L 2 (Ba) (again invoking (48)). Choosing a 0 such that B a Ω for any a ≤ a 0 , we therefore have for any a ≤ a 0 .
Regarding item (ii), the integral operator L c a (given by (26)) 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. [22,Thm. 6.1.12]. Moreover, these kernels being O(1) as a → 0 (see the kernel expansions given after (40) 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 (48), we have 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 (16)  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: Concluding, I − L a is invertible and its inverse, given by ( Proof of Lemma 7. We first derive an explicit expression of γ a , given by (68), recalling the definition (36) of U a . First, using the rescaled form (25) 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 (39a)+a(39b)+a 2 (44a)+a 3 (44b)+a 4 (44c) of the VIEs governing the coefficients U 0 , . . . , U 4 , I − H)U a ], using the resulting formula for evaluating P a γ a with γ a as given by (68), and rearranging terms, we obtain Using estimates (70a), (70b) and (70c) in (69) therefore yields P a γ a H 1 (B) ≤ a 5 C for some C > 0. Lemma 7 finally follows from applying relations (48) to the latter estimate of P a γ a .
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 (33), as shown in Appendix C.
It is then a simple matter to check, using (71), that the solution u B = U 1 ] of the FSTP for the polynomial background field u(x) = E 1 ·x + E 2 : (x⊗x) verifies formulas (32).
Second identity: by Lemma 2 and (46b) by interrelations (73a,d) by VIE (31) = k 2 g z ·(B 1p − Q 1p )·E p Third identity: We have where the second equality follows from identities 73) and the third uses definitions (47a) of B 21 and (47b) of Q 21 . Now, expressing X B by means of equation (46a) and using the interrelation (73a), we have so that the last three terms inside the square brackets of (74) are given by Using the above equality in (74) finally yields the desired identity.