Bojan B Guzina
email: guzina@wave.ce.umn.edu

Marc Bonnet
email: bonnet@lms.polytechnique.fr

Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics

Keywords: Topological derivative, Helmholtz equation, inverse scattering, transmission problem, acoustic waves, penetrable obstacles, probe method

The aim of this study is an extension and employment of the concept of topological derivative as it pertains to the nucleation of infinitesimal inclusions in a reference (i.e. background) acoustic medium. The developments are motivated by the need to develop a preliminary indicator functional that would aid the solution of inverse scattering problems in terms of a rational initial "guess" about the geometry and material characteristics of a hidden (finite) obstacle; an information that is often required by iterative minimization algorithms. To this end the customary definition of topological derivative, that quantifies the sensitivity of a given cost functional with respect to the creation of an infinitesimal hole, is adapted to permit the nucleation of a dissimilar acoustic medium. On employing the Green's function for the background domain, computation of topological sensitivity for the three-dimensional Helmholtz equation is reduced to the solution of a reference, Laplace transmission problem. Explicit formulas are given for the nucleating inclusions of spherical and ellipsoidal shape. For generality the developments are also presented in an alternative, adjoint-field setting that permits nucleation of inclusions in an infinite, semi-infinite or finite background medium.

Through numerical examples it is shown that the featured topological sensitivity could be used, in the context of inverse scattering, as an effective obstacle indicator through an assembly of sampling points where it attains pronounced negative values. On varying a material characteristic (density) of the nucleating obstacle, it is also shown that the proposed methodology can be used as a preparatory tool for both geometric and material identification.

Introduction

Shape reconstruction of obstacles embedded in an acoustic medium, using either faror near-field patterns of the scattered wave field, is a challenging subject [START_REF] Colton | Inverse acoustic and electromagnetic scattering teory[END_REF][START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF] with application to diverse areas such as sonar detection and medical imaging. In view of their non-linear nature, inverse problems of this class are often dealt with by way of gradient-based optimization [START_REF] Colton | Numerical solution of an inverse scattering problem for acoustic waves[END_REF][START_REF] Farhat | On the solution of three-dimensional inverse obstacle acoustic scattering problems by a regularized newton method[END_REF], sometimes aided by the adjoint-field sensitivity estimates [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] for computational efficiency and level-set methodology [START_REF] Litman | Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set[END_REF] as a means to permit changes in topology of the trial scatterer. Unfortunately, cost functionals associated with shape reconstruction problems are often non-convex [START_REF] Pike | Scattering: scattering and inverse scattering in pure and applied science[END_REF] which renders the success of gradient-based optimization highly dependent on the initial "guess" in terms of the location, topology, and geometry of the scatterer. In the absence of suitable prior information, a direct way to deal with the problem revolves around the use of global search algorithms [START_REF] Michalewicz | How to solve it: modern heuristics[END_REF][START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]. Owing to their exceeding computational cost, however, such solutions are by rule impractical which exposes the need for computationally effective, preliminary global-search tools whose primary purpose would be to generate a reliable initial "guess" for gradient-based optimization.

Building on the results in shape optimization obtained for Laplace [START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF][START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF] and Helmholtz [START_REF] Samet | The topological asymptotic for the Helmholtz equation[END_REF][START_REF] Pommier | The topological asymptotic for the Helmholtz equation with Dirichlet condition on the boundary of an arbitrarily shaped hole[END_REF] systems, [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF][START_REF] Bonnet | Topological derivative for the inverse scattering of elastic waves[END_REF][START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF] and [START_REF] Feijoo | A new method in inverse scattering based on the topological derivative[END_REF] have recently established the method of topological sensitivity as a tool for preliminary, grid-based reconstruction of obstacles in the context of inverse elastic and acoustic scattering. In the approach the topological derivative, which quantifies the sensitivity of a given cost functional with respect to the nucleation of an infinitesimal obstacle in the reference (background) medium, is used as a spatial obstacle indicator. Notwithstanding their usefulness, however, the foregoing topological sensitivity analyses are limited in the sense that they are focused on the nucleation of impenetrable scatterers with either Dirichlet or Neumann conditions imposed on their boundary. To deal with elastic-wave identification of penetrable obstacles, in terms of their both geometric and material features, [START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF] have recently developed a generalization of topological sensitivity that postulates the nucleation of dissimilar elastic inclusions. In this study, the latter idea is carried out to deal with inverse scattering problems in acoustics. In particular, the proposed "material-topological" sensitivity is shown to consist of a monopole term, related to the compressibility contrast, and a dipole term involving the mass density contrast. For generality, the proposed developments are also cast within the adjoint-field formulation, an alternative framework allowing efficient computation of the material-topological sensitivity in an arbitrary (infinite or finite, homogeneous or heterogeneous) background acoustic medium. Through numerical examples it is shown that the material-topological sensitivity can be used, in the context of inverse scattering, as an effective obstacle indicator through an assembly of sampling points where it attains pronounced negative values. On varying the material characteristics of a nucleating obstacle, the proposed indicator is also shown to permit both preliminary geometric and material identification. The latter result may be for instance useful in breast cancer detection wherein the mechanical characteristics of a lesion, detected e.g. via ultrasound or magnetic resonance imaging, may allow one to differentiate between malignant and benign growths [START_REF] Sarvazyan | Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics[END_REF][START_REF] Fatemi | Probing the dynamics of tissue at low frequencies with the radiation force of ultrasound[END_REF].

For generality, it is noted that the idea of (preliminary) obstacle reconstruction via spatial sampling of a given indicator functional is shared by a diverse array of inverse scattering techniques, now commonly referred to as the "sampling" or "probe" methods [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF], which notably include the linear sampling method, e.g. [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF][START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF]. In this context, an interesting observation is that the linear sampling and topological sensitivity methods both revolve around the evaluation of a probing functional that projects the relevant (background) Green's function onto a suitable function space synthesizing the experimental observations. However, these two methods have distinct theoretical frameworks, and their relationship is very much an open research issue. As one of the key distinctions relevant to this study, the derivation of materialtopological sensitivity formulae for a given cost functional requires an asymptotic expansion of the scattered field induced by the nucleating penetrable obstacle. Although the related asymptotic expressions have been proposed elsewhere (e.g. [START_REF] Ammari | Reconstruction of the small inhomogeneities from boundary measurements[END_REF]) for inverse electromagnetic problems and penetrable scatterers, they have so far been utilized for the explicit reconstruction of small inhomogeneities based on either i) treatment of the measurement residuals [START_REF] Ammari | Electromagnetic scattering by small dielectric inhomogeneities[END_REF][START_REF] Vogelius | Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter[END_REF][START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF], or ii) the reciprocity-gap approach [START_REF] Ammari | Boundary integral formulae for the reconstruction of electric and electromagnetic inhomogeneities of small volume[END_REF]. In contrast, the present approach, that makes use of the spatial distribution of topological sensitivity as an indicator functional, aims at preliminary reconstruction of finite obstacles and may broadly be categorized as a sampling technique exploiting asymptotic expansions.

Preliminaries

With reference to the Cartesian frame {O; ξ 1 , ξ 2 , ξ 3 }, consider a semi-infinite acoustic domain Ω = {ξ | ξ 3 < 0} housing a penetrable obstacle B that is bounded by the closed, piecewise-smooth surface S (Fig. 1). The reference background medium and the obstacle are each assumed to be homogeneous with wave speed and mass density (c, ρ) and (c * , ρ * ), respectively. In this setting, Ω -= Ω\(B ∪ S) is the exterior region surrounding the obstacle, and Ω-= Ω -∪ Σ ∪ S denotes the closure of Ω -where Σ = {ξ | ξ 3 = 0} is the "top" surface of the half-space. For further reference the normal on Σ ∪ S, oriented outward from Ω -, will be denoted by n. The homogeneous boundary conditions on Σ are assumed to be either of the Dirichlet or Neumann type. With the implicit timeharmonic factor e iωt omitted henceforth for brevity, let the obstacle be illuminated by the free (i.e. incident) acoustic field p F defined as the response of an obstacle-free half-space due to prescribed source distribution f Ω (ξ) and boundary excitation f Σ (ξ) so that

∇ 2 p F + k 2 p F + f Ω = 0, ξ ∈ Ω, α p F + (1-α) p F ,n = f Σ , ξ ∈ Σ. (1)
Here k = ω/c; g ,n = n • ∇g, and α takes the respective values 0 and 1 when Σ is the Neumann and Dirichlet boundary. In what follows, it is assumed that ω is not a fictitious eigenfrequency of any of the featured boundary-value problems. When an obstacle is present, the prescribed excitation f ≡ (f Ω , f Σ ) gives rise to the acoustic field p in Ω -∪ B that can be conveniently decomposed as

p(ξ) = p F (ξ) + p(ξ), ξ ∈ Ω - (2)
in the exterior domain where p denotes the perturbation component, termed the scattered field. Assuming that the free field p F is known beforehand, the forward scattering problem, i.e. the task of calculating p in Ω -∪ B (given f ) can be cast in terms of p| Ω -and p| B as

∇ 2 p + k 2 p = 0, ξ ∈ Ω -, ∇ 2 p + γ 2 k 2 p = 0, ξ ∈ B, α p + (1-α) p,n = 0, ξ ∈ Σ p + p F = p, p,n + p F ,n = β p ,n ξ ∈ S, (3)
where β = ρ/ρ * , γ = c/c * , and the normal derivatives p,n , p F ,n are defined in terms of the same normal n. To ensure physical relevance, the scattered field is further assumed to satisfy the regularity and radiation conditions

p = O 1 R , p,R + ik p = o 1 R , R ≡ |ξ| → ∞, ξ ∈ Ω -. (4)
For identification purposes, let Γ obs ⊂ Σ ∪ Ω -denote the measurement surface. With the foregoing definitions, the inverse problem of interest can be stated as a task of resolving the "true" obstacle B true , both in terms of its geometry and material characteristics (c * ,true , ρ * ,true ), from the knowledge of: i) time-harmonic source f used to illuminate the obstacle, ii) acoustic properties (c, ρ) of the background medium, and iii) distribution of the induced pressure field, p true , as observed over the measurement surface Γ obs . In what follows, these measurements will be denoted by p obs , so that p obs (ξ) = p true (ξ), ξ ∈ Γ obs under ideal modeling and measurement conditions. For generality, it is assumed that either Γ obs ∩ Σ ≡ Γ obs Σ or Γ obs ∩ Ω -≡ Γ obs Ω (see Fig. 1) may be an empty set.

For a systematic treatment of the identification problem, a cost functional is established that quantifies the misfit between experimental observations p obs and their acoustic predictions p, calculated for a trial obstacle B. To provide a focus in the study, the ensuing discussion deals with the class of cost functionals

J (Ω -, β, γ; f ) = Γ obs ϕ (p(ξ), ξ) dΓ ξ , (5)
where β and γ synthesize the material characteristics of a trial obstacle, while the distance function ϕ is assumed to be real-valued, nonnegative, and differentiable with respect to the real and imaginary parts of its first argument. In waveform tomography, ϕ commonly takes the weighted least-squares format

ϕ(p(ξ), ξ) = 1 2 W (ξ) {p(ξ)-p obs (ξ)} p(ξ)-p obs (ξ) , W (ξ) > 0. (6)

Green's functions

With an implicit time-harmonic factor e iωt as examined earlier, let

G(ξ, x; k) = e -ikr 4πr , H(ξ, x; k) ≡ n•∇G(ξ, x; k) = - e -ikr 4πr 2 (1+ikr) r ,n , (7)
where r = |ξ-x| and ∇G signifies the gradient of G with respect to its first argument, denote the acoustic fundamental solution for the free space with wavenumber k, so that

∇ 2 G(ξ, x; k) + k 2 G(ξ, x; k) + δ(ξ -x) = 0, ξ ∈ R 3 . ( 8
)
On the basis of ( 7), the Green's function for the reference (i.e. obstacle-free) half-space Ω can be written as

Ĝ(ξ, x; k) = G(ξ, x; k) + (1-2α) G(ξ, x ′ ; k), Ĥ(ξ, x; k) = H(ξ, x; k) + (2α-1)H(ξ, x ′ ; k), (9)
where x and x ′ are symmetric with respect to the top surface of the half-space, and α ∈ {0, 1} follows the convention introduced in (1). By virtue of the foregoing definitions, the free field satisfying (1) can be written as

p F (x) = Ω - f Ω (ξ) Ĝ(ξ, x, k) dV ξ + Σ f Σ (ξ) α Ĥ(ξ, x, k)+(1-α) Ĝ(ξ, x, k) dΓ ξ . ( 10
)

Generalized Topological Sensitivity

This section is focused on the development of topological derivative for the class of cost functionals J given by [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF]. To this end, let B ⊂ R 3 be a fixed bounded open set with volume |B| containing the origin. For generality, it is assumed that S = ∂B is piecewisesmooth with components of class C 2 . With such definitions, one may consider the introduction of a small penetrable obstacle B ǫ (x o ) = x o + ǫB with material properties c * = c/γ, ρ * = ρ/β and size ǫ > 0, at a fixed trial point x o in an otherwise obstacle-free half-space Ω with properties c and ρ. Following [START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF][START_REF] Céa | The shape and topological optimization connection[END_REF], one is in particular interested in the asymptotic behavior of J (Ω - ǫ , β, γ; f ) for infinitesimal ǫ > 0, where Ω - ǫ = Ω\ Bǫ (x o ) and Bǫ = B ǫ ∪S ǫ is the closure of B ǫ (x o ) ⊂ Ω. With reference to this limiting behavior, the topological derivative T (x o , β, γ; f ) of the cost functional [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] for an obstacle-free body can be defined through the expansion

J (Ω - ǫ , β, γ; f ) = J (Ω; f ) + T (x o , β, γ; f ) h(ǫ) + o(h(ǫ)), ǫ → 0, (11)
where J (Ω; f ) ≡ J (Ω, β, γ; f ) denotes the degenerate value of (5) when p = p F and h(ǫ) > 0, characterizing the leading term, is to be determined. For the concept of topological derivative to make sense, it is assumed that

lim ǫ→0 h(ǫ) = 0, |T (x o , β, γ; f )| < ∞, x o ∈ Ω. ( 12
)
In general, the functional form of h(ǫ) depends on the nature of the governing field equation, the type of boundary conditions prescribed on S , and dimensionality of the problem [START_REF] Céa | The shape and topological optimization connection[END_REF][START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF]. As examined earlier, the focus of this study is the expansion [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF], and in particular its leading-term coefficient T , for the three-dimensional Helmholtz equation in the context of cost functionals (5) and a penetrable obstacle problem. With reference to (2) and ( 5), evaluation of J (Ω - ǫ , β, γ; f ) requires the knowledge of the acoustic solution, p = pǫ , to the forward transmission problem (3) wherein B is replaced by B ǫ ≡ B ǫ (x o ). Since the scattered field pǫ is expected to vanish as ǫ → 0 whereas the free field p F does not depend on ǫ, the leading contribution to J (Ω - ǫ , β, γ; f ) can be sought on the basis of an expansion with respect to pǫ , i.e.

J (Ω - ǫ , β, γ; f ) = J (Ω; f ) + Γ obs Re ∂ϕ ∂p p F (ξ), ξ pǫ (ξ) dΓ ξ + O( pǫ 2 ), ( 13
)
where

∂ϕ ∂p ≡ ∂ϕ ∂p R -i ∂ϕ ∂p I p R = Re(p) , p I = Im(p) .
By means of ( 11) and ( 13), the topological derivative of J can be recast as

T (x o , β, γ; f ) = lim ǫ→0 1 h(ǫ) Γ obs Re ∂ϕ ∂p p F (ξ), ξ pǫ (ξ) dΓ ξ . (14)
The direct approach for the computation of topological derivative proposed in this study represents a generalization of the methodology developed in [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] for the inverse scattering problems in infinite and semi-infinite media. It entails a substitution of the leading contribution of pǫ into (14) and seeking the limit of the resulting expression as ǫ → 0. To this end, it is useful to recall the integral representation

pǫ (x) = Sǫ pǫ ,n (ξ) Ĝ(ξ, x, k) dS ξ - Sǫ pǫ (ξ) Ĥ(ξ, x, k) dS ξ , x ∈ Ω - ǫ ( 15
)
of the scattered field induced by B ǫ . On employing the interfacial conditions over S = S ǫ as in (3), the divergence theorem, and the Taylor expansion of the featured Green's function at ξ = x o , the limiting behavior of (15) for a vanishing obstacle can be written as

pǫ (x) = (1-β) Bǫ ∇p ǫ (ξ) dV ξ •∇ Ĝ(x o , x, k) -(1-βγ 2 ) k 2 Bǫ p ǫ (ξ) dV ξ Ĝ(x o , x, k) + o(ǫ 3 ), x ∈ Ω - ǫ , ǫ → 0 (16)
One may note that (16) requires the knowledge of p ǫ and ∇p ǫ in the interior of a vanishing obstacle, quantities that are unknown beforehand. Their evaluation is addressed next.

Small-inclusion asymptotics

The forward problem (3) for B = B ǫ can be conveniently formulated in terms of a pair of boundary integral equations

κ(x) pǫ (x) + Sǫ pǫ (ξ) Ĥ(ξ, x, k) dS ξ - Sǫ pǫ ,n (ξ) Ĝ(ξ, x, k) dS ξ = 0, κ(x) pǫ (x) - Sǫ pǫ (ξ)H(ξ, x, γk) dS ξ + 1 β Sǫ pǫ ,n (ξ)G(ξ, x, γk) dS ξ = -κ(x) p F (x) + Sǫ p F (ξ)H(ξ, x, γk) dS ξ - 1 β Sǫ p F ,n (ξ)G(ξ, x, γk) dS ξ ,                  x ∈ S ǫ , (17)
written respectively for the exterior and interior region in terms of the scattered field over S ǫ . Here κ(x) is the free term, equal to 1 2 if B ǫ is smooth at x, and n is again oriented outward from Ω -. One may note that integral equations ( 17) are free of Cauchy principal values owing to the particular structure of the acoustic fundamental solution and the earlier assumption that S ǫ is piecewise-smooth with components of class C 2 [START_REF] Bonnet | Boundary integral equations methods for solids and fluids[END_REF].

In this setting, the limiting form of p ǫ and ∇p ǫ in B ǫ as ǫ → 0 is sought by performing an asymptotic analysis of [START_REF] Feijoo | A new method in inverse scattering based on the topological derivative[END_REF], following e.g. the procedures described in [START_REF] Ammari | Electromagnetic scattering by small dielectric inhomogeneities[END_REF], and the invoking interfacial conditions (3d), to obtain

p ǫ (ξ) = p F (x o ) + o(1), ∇p ǫ (ξ) = ∇p F (x o )•∇ ζ g(ζ) + o(1), ξ ∈ B ǫ , ζ ∈ B. (18)
Here ζ = (ξx o )/ǫ is the scaled position vector, and g(ζ) is a component of the vector solution pair (g, g) that satisfies the (vector) Laplace transmission problem for the "unit" obstacle B = B ǫ | ǫ=1 embedded in an infinite medium:

∇ 2 ζ g = 0 (ζ ∈ B), ∇ 2 ζ g = 0 (ζ ∈ R 3 \ B), g + ζ = g, g,η + η = βg ,η , (ζ ∈ S ). (19)
In [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF], η(ζ) = n(ξ) is the unit normal on S , oriented toward the interior of B; the normal derivatives g,η and g ,η are both defined in terms of η, and the exterior field g decays as O(|ζ| -1 ) or faster at infinity.

On substituting ( 18) into ( 16), one finds that the scattered field for a vanishing inclusion can be expressed as

pǫ (x) = ǫ 3 |B| (1-β)∇p F (x o )•A•∇ Ĝ(x o , x, k) -(1-βγ 2 ) k 2 p F (x o ) Ĝ(x o , x, k) + o(ǫ 3 ), x ∈ Ω - ǫ , ǫ → 0 ( 20
)
where A is the (constant, second-order) polarization tensor

A = |B| -1 B ∇ ζ g dV ζ = -|B| -1 S g ⊗ η dS ζ , (21)
whose latter expression stems from the divergence theorem. On the basis of ( 19) and ( 21), it is noted that the polarization tensor depends only on the shape B and relative density β of the nucleating obstacle.

One may also observe from (20) that pǫ (x) behaves as O(ǫ 3 ) in Ω -as ǫ → 0 and thus, by way of ( 14), h(ǫ) ∝ ǫ 3 regardless of the shape of B ǫ . Since the choice of the multiplicative constant is, through [START_REF] Colton | Inverse acoustic and electromagnetic scattering teory[END_REF], arbitrary, it is assumed in this study that

h(ǫ) = ǫ 3 |B|, (22)
i.e. that h(ǫ) is given by the measure of a vanishing obstacle domain. On the basis of ( 20) and ( 22), the formula for topological derivative ( 14) reduces to

T (x o , β, γ; f ) = Γ obs Re ∂ϕ ∂p p F (ξ), ξ (1-β) ∇p F (x o )•A•∇ Ĝ(x o , ξ, k) -(1-βγ 2 ) k 2 p F (x o ) Ĝ(x o , ξ, k) dS ξ , x o ∈ Ω, ( 23
)
where, for the least-squares-type cost functional [START_REF] Bonnet | Boundary integral equations methods for solids and fluids[END_REF],

∂ϕ ∂p p F (ξ), ξ = W (ξ) (p F (ξ)-p obs (ξ)).
Here it is noted that ( 23) consists of a dipole and a monopole term, which vanish respectively when the mass density and compressibility of the obstacle equals that of the background medium. The limiting case of an acoustically-hard vanishing obstacle (exterior Neumann problem) is obtained by taking the limit ρ * → ∞, i.e. by setting β = 0 in [START_REF] Hobson | The theory of spherical and ellipsoidal harmonics[END_REF] (see also [START_REF] Dassios | Low frequency scattering[END_REF] in the context of general acoustics). In contrast, the case of an acoustically-soft vanishing obstacle (β → ∞) leads to an asymptotic behavior that is distinct from [START_REF] Hobson | The theory of spherical and ellipsoidal harmonics[END_REF] and hence not recoverable by the present approach. The key reason for such distinction is that the asymptotic analysis underpinning [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] and ( 22) is not appropriate when the Dirichlet condition is assumed along the boundary of a vanishing obstacle.

A generalization of ( 5) and ( 23) to situations involving Q sequentially-applied acoustic sources f q is straightforward and involves external summation of the form

J f (•) ≡ Q q=1 J (•; f q ) and T f (•) ≡ Q q=1 T (•; f q ). ( 24
)

Explicit Expressions

In what follows, explicit expressions for generalized topological sensitivity are given for the canonical cases of spherical and ellipsoidal inclusions. By their nature, the ensuing results are closely related to the low-frequency approximation of the acoustic scattering problem, a subject dating back to [START_REF] Rayleigh | On the incidence of aerial and electric waves upon small obstacles in the form of ellipsoids or elliptic cylinders and on the passage of electric waves through a circular aperture in a conducting screen[END_REF] and more recently covered in [START_REF] Dassios | Low frequency scattering[END_REF].

Spherical obstacle

The simplest specialization of the foregoing developments is the case when the vanishing obstacle is spherical, i.e. when S is the unit sphere and |B| = 4π/3. For this geometry, ( 19) can be easily solved using spherical harmonics, to obtain

g(ζ) = 3 β + 2 ζ, ζ ∈ B. (25)
The polarization tensor ( 21) is accordingly found, on the basis of ( 25), as

A sph = 3 β + 2 I 2 , (26)
where I 2 is the second-order identity tensor. It is interesting to note that in dealing with 3D electromagnetic scattering, [START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF] and [START_REF] Ammari | Recovery of small inhomogeneities from the scattering amplitude at a fixed frequency[END_REF] list apparently inconsistent expressions for the polarization tensor, expressible respectively as

A sph = 2 β -1 2β + 1 I 2 , A sph = 8π 1 β + 1 I 2 , (27)
using present notation. On comparing [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF] in this paper with the last formula on page 195 in [START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF] and the last formula on page 883 in [START_REF] Ammari | Recovery of small inhomogeneities from the scattering amplitude at a fixed frequency[END_REF], it can be shown that β in [START_REF] Pommier | The topological asymptotic for the Helmholtz equation with Dirichlet condition on the boundary of an arbitrarily shaped hole[END_REF] takes the meaning of either relative magnetic permeability µ ⋆ /µ or relative electric permittivity ǫ ⋆ /ǫ.

Ellipsoidal obstacle To examine the behavior of ( 23) when the vanishing obstacle is non-spherical, one may consider the solution to the Laplace transmission problem [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF] in situations when S is an ellipsoid with principal directions e ′ 1 , e ′ 2 , e ′ 3 and semi-axes a 1 , a 2 , a 3 . Here the characteristic volume is |B| = 4πa 1 a 2 a 3 /3 which is inherently accounted for through the definition of h(ǫ). In this setting, the analytical solution of ( 19) is sought by introducing an ellipsoidal coordinate system and expanding g, g in terms of spheroidal harmonics [START_REF] Hobson | The theory of spherical and ellipsoidal harmonics[END_REF]. This procedure yields the solution for g as

g(ζ) = 3 m=1 1 1-(1-β) I m (e ′ m ⊗e ′ m )•ζ, ζ ∈ B ( 28
)
where I m are the functions of the two aspect ratios of B (e.g. a 2 /a 1 and a 3 /a 1 ) given by

I m = a 1 a 2 a 3 2 ∞ 0 dr (r + a 2 m ) r+a 2 1 r+a 2 2 r+a 2 3 , m = 1, 2, 3
On the basis of ( 18) and ( 28), the polarization tensor featured in [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] takes the form

A ell = 3 m=1 e ′ m ⊗e ′ m 1 -(1-β)I m ( 29
)
that is diagonal in the principal axes frame (e ′ 1 , e ′ 2 , e ′ 3 ), albeit not isotropic as in [START_REF] Pike | Scattering: scattering and inverse scattering in pure and applied science[END_REF]. To demonstrate that A ell remains finite for arbitrarily chosen relative inclusion density β ≥ 0 and semi-axes a 1 , a 2 , a 3 , it can be shown by elementary calculations that I (1) +I (2) + I (3) = 1 which, through restriction I m > 0, yields I m < 1 and thus the desired result. Moreover, it is worth noting that on setting a 1 = a 2 = a 3 = 1 in (28), one obtains I m = 1/3 (m = 1, 2, 3) which directly reduces the polarization tensor [START_REF] Rayleigh | On the incidence of aerial and electric waves upon small obstacles in the form of ellipsoids or elliptic cylinders and on the passage of electric waves through a circular aperture in a conducting screen[END_REF] to its spherical counterpart [START_REF] Pike | Scattering: scattering and inverse scattering in pure and applied science[END_REF]. Depending on the quantity of interest, however, such a reduction in the context of acoustic scattering may not always involve a trivial calculation [START_REF] Dassios | Low frequency scattering[END_REF].

Qualitative obstacle identification via generalized topological sensitivity

Equation ( 23) provides, when applied to the featured cost functional J, an information indicating whether J increases (T (x o ) > 0) or decreases (T (x o ) < 0) in response to the nucleation at x o of an infinitely small penetrable scatterer with prescribed shape and material characteristics. Thus, the topological sensitivity field T (x o ) is guaranteed to define a correct obstacle indicator only in the limit as the size of a hidden obstacle approaches zero. It is nonetheless natural to examine whether the spatial distribution of T can still furnish a useful information for the reconstruction and characterization of finite obstacles. Motivated by the ultimate goal of minimizing J , the key idea in this setting is to approximate the support of a hidden (finite) obstacle via regions where T attains pronounced negative values. While the reasoning behind such an idea is heuristic as it lacks a rigorous mathematical backing such as that underpinning the linear sampling method [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF][START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF], the numerical experiments presented in this study as well as the other (impenetrable obstacle) results in earlier works [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF][START_REF] Bonnet | Topological derivative for the inverse scattering of elastic waves[END_REF][START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF][START_REF] Feijoo | A new method in inverse scattering based on the topological derivative[END_REF] support the utility of the proposed notion as a preliminary reconstruction tool. Since the computation of the indicator function T (over the volume of interest) is significantly faster than an iterative (e.g. minimization-based) inversion, this approach offers a convenient means of extracting preliminary obstacle information from the data p obs . Such information can then be used either in a stand-alone manner for limited-accuracy reconstruction, or as a reliable initial "guess" for more elaborate, iterative reconstruction schemes.

Discussion

For a broader perspective of the foregoing results, this section deals with complementary developments that include a discussion of the polarization tensor, an adjoint-field approach to topological sensitivity, and a treatment of finite background domains.

Nature of the polarization tensor

With the help of the Green's first identity applied to harmonic fields gk in R 3 \ B, g k in B, and ζ k in B (k = 1, 2, 3), where g = g k e k and g = gk e k are the respective interior and exterior parts of the solution to the Laplace problem [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF], one obtains the identity

R 3 \B gi,k gj,k dV ζ + B βg i,k g j,k + δ ik dV ζ = S gi gj,η -βg i g j,η -ζ i η j dS ζ . (30)
On the basis of the interfacial conditions in [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF], ( 21), [START_REF] Samet | The topological asymptotic for the Helmholtz equation[END_REF], and the Green's second identity

S ζ i g j,η dS ζ = S η i g j dS ζ ,
one finds that the Cartesian components of A are such that

A ij + βA ji = R 3 \B gi,k gj,k dV ζ + B βg i,k g j,k + δ ik dV ζ . (31)
The right-hand side of (31) clearly defines a symmetric and positive-definite secondorder tensor. Hence one must have A ij + βA ji = A ji + βA ij , i.e. A ij = A ji (β = 1). Equation ( 31) thus shows that (1 + β)A, and hence the polarization tensor A, is symmetric and positive-definite regardless of the shape of B. This is a useful result as it indicates that the leading-term scattering effects [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], and thus the topological sensitivity [START_REF] Hobson | The theory of spherical and ellipsoidal harmonics[END_REF], for a connected infinitesimal obstacle of arbitrary shape might be cast in terms of its ellipsoidal counterpart. Whether such a possibility can be materialized, however, depends on the solvability of the non-linear system

(1 -β) I m a 2 a 1 , a 3 a 1 = 1 -λ -1 m , m = 1, 2, 3 (32)
in terms of parameters a 2 /a 1 , a 3 /a 1 and β of an equivalent-ellipsoidal vanishing obstacle where λ m are the (real-valued) eigenvalues of A.

Direct versus adjoint field formulation

For generality it is instructive to recast [START_REF] Hobson | The theory of spherical and ellipsoidal harmonics[END_REF], obtained via the so-called direct approach, within the framework of the adjoint field method that is commonly used in the shape sensitivity (e.g. [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF][START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF]) and topological sensitivity (e.g. [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF]) analyses. In the present setting, the adjoint field is defined as an auxiliary pressure field

⋆ p(x) = Γ obs ∂ϕ ∂p (p F (ξ), ξ) Ĝ(x, ξ, k) dS ξ , (33)
signifying the response of the obstacle-free domain Ω due to virtual excitation

g(ξ) = ∂ϕ ∂p (p F (ξ), ξ), ξ ∈ Γ obs , (34)
written in terms of the misfit function ϕ for p = p F . On the basis of ( 14), ( 20), [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF] and [START_REF] Vogelius | Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter[END_REF], an adjoint-field variant of formula [START_REF] Hobson | The theory of spherical and ellipsoidal harmonics[END_REF] for material-topological sensitivity is obtained as

T (x o , β, γ; f ) = Re (1-β) ∇p F •A•∇ ⋆ p -(1-βγ 2 ) k 2 p F ⋆ p (x o ), x o ∈ Ω. ( 35
)
One may note that the difference between ( 23) and ( 35) is at this point strictly formal.

In particular, for all background domains Ω for which the Green's function is available, expressions ( 23) and ( 35) lead to the same Green's function-based computational scheme for the evaluation of T . For more complex forms of Ω (e.g. finite domains of arbitrary shape), on the other hand, practical unavailability of suitable Green's function requires a treatment wherein the direct and adjoint field formulations lead to distinct computational procedures. This issue is addressed next.

Finite background domains

To deal with acoustic obstacle reconstruction in the context of finite bodies, consider a background domain Ω with properties (c, ρ) and external boundary Σ = ∂Ω, housing a penetrable obstacle B with properties (c * , ρ * ). Following the previously established convention, Ω -= Ω\B, S = ∂B, Ω-= Ω -∪ Σ ∪S, while n and n ′ denote the respective normals on Σ ∪ S (oriented outward from Ω -) and Γ obs Ω . On writing Σ = Σ D ∪ Σ N , it is assumed that Ω -is subject to a domain source distribution and to Dirichlet and Neumann boundary data prescribed respectively over Σ D and Σ N . Next, let the Green's function Ǧ solve [START_REF] Bonnet | Topological derivative for the inverse scattering of elastic waves[END_REF] over Ω subject to homogeneous boundary conditions

Ǧ(ξ, x, k) = 0 (ξ ∈ Σ D ), Ȟ(ξ, x, k) = 0 (ξ ∈ Σ N ).
On superseding Ĝ and Ĥ respectively by Ǧ and Ȟ, the scattered field pǫ is, for a finite body Ω - ǫ , still given by the integral representation ( 16) and the asymptotic expansion [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] still holds [START_REF] Ammari | Electromagnetic scattering by small dielectric inhomogeneities[END_REF] with the polarization tensor A again given by ( 21) and [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF]. As a result, expression [START_REF] Hobson | The theory of spherical and ellipsoidal harmonics[END_REF] for generalized topological sensitivity (where Ĝ is replaced with Ǧ) and its adjoint-field variant [START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF] still hold.

The analytical expression of the latter Green's function is, however, not available except for very few geometrically-simple domains Ω. Assuming that the polarization tensor A and the free field p F are calculated beforehand, the computation of topological derivative for a finite body using the direct approach ( 23) consequently entails solving one boundary value problem per sampling point x o (as a means to evaluate the weighted sum of Ǧ(x o , ξ, k) and ∇ Ǧ(x o , ξ, k)).

This potentially demanding computational procedure can be significantly reduced by resorting to the adjoint field approach. Indeed, upon invoking the definition of Ǧ(ξ, x, k) and the symmetry property Ǧ(x, ξ, k) = Ǧ(ξ, x, k), one can show that the adjoint field ⋆ p(x), given by [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF], solves the transmission problem

∇ 2 ⋆ p + k 2 ⋆ p = 0 (ξ ∈ Ω\Γ obs Ω ), [[ ⋆ p]] = 0, [[ ⋆ p ,n ′ ]] = g (ξ ∈ Γ obs Ω ), ⋆ p = 0 (ξ ∈ Σ D ), ⋆ p ,n = 0 (ξ ∈ Σ N \Γ obs Σ ), ⋆ p ,n ′ = g (ξ ∈ Γ obs Σ ), (36)
where the single-layer excitation g is given by [START_REF] Vogelius | Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter[END_REF], n ′ is the unit normal on Γ obs Ω , and

[[w]] = lim τ →0 w(ξ+τ n ′ )-w(ξ-τ n ′ ), ξ ∈ Γ obs
Ω denotes the jump in w across Γ obs Ω . Hence, the adjoint field approach (35) requires the solution of only one additional boundary value problem (used to compute ⋆ p and ∇ ⋆ p) that is common for all sampling points. In this sense, the adjoint field methodology represents the most effective (if not direct) route for computing the topological sensitivity in a general background-body configuration. As such, it will be used as a computational basis for the ensuing examples.

Results

To examine the effectiveness of generalized topological sensitivity as a tool for preliminary obstacle reconstruction and identification, a set of numerical results is presented next. In the sequel, each synthetic testing configuration consists of M measurement points x m (m = 1, . . . , M ). For reconstruction purposes, the obstacle (or a set thereof) is illuminated by a set of Q point sources of magnitude ρc 2 ℓ, applied in sequence at locations x q (q = 1, . . . , Q) where ℓ is the reference length. On the basis of ( 5), ( 6) and ( 24), the generalized topological sensitivity is thus computed with reference to the least-squares cost function

J f (Ω -, β, γ) = 1 2 Q q=1 M m=1 p q (x m ) -p obs q (x m ) p q (x m ) -p obs q (x m ) , (37)
where p obs q and p q denote respectively the measured (i.e. "true") and trial acoustic fields induced at sensor locations x m by the source acting at x q .

In what follows, the synthetic data p obs q in (37) are generated for five "true" obstacle configurations: i) one penetrable obstacle (β true = 2, γ true = 0.5); ii) one rigid obstacle (β true = 0, γ true = 1); iii) two penetrable obstacles; iv) two rigid obstacles, and v) two dissimilar obstacles, one penetrable and one rigid. For brevity, these configurations are labeled as P, R, PP, RR and PR, respectively. All obstacles are ellipsoidal, their semiaxes being aligned with the reference Cartesian frame as shown in Fig. 2. The geometric and material parameters for all five obstacle configurations are listed in Table 1. With the exception of three-dimensional plots, all ensuing T f -distributions are plotted with reference to the "vertical" plane ξ 2 = 0.5ℓ. Accordingly, the horizontal and vertical axes of sectional distributions in Figs. 3-8 and 10-14 carry implicit labels ξ 1 /ℓ and ξ 3 /ℓ, respectively. To provide a basis for comparison, intersection of the featured plane (ξ 2 = 0.5ℓ) with "true" scatterers is indicated in white in the figures. The scatterers are embedded in either acoustic half-space (H) or the full-space (F). Building upon the previously adopted notation, the scatterer-background medium combinations are designated by a hyphenated symbol, with e.g. P-F denoting a single penetrable obstacle in the acoustic full-space. The synthetic data p obs q are computed by means of a direct boundary element method (BEM), wherein each ellipsoidal obstacle is meshed using 384 eight-node elements. In this way, the BEM meshes feature at least 15 nodes per wavelength for the highest wavenumber considered, kℓ = 8.

The "complete" (i.e. the most extensive) testing configuration used in the examples consists of an array of points located at the vertices of six uniform 10×10 rectangular grids placed along the faces of a box defined by the planes ξ 1 = ±6ℓ, ξ 2 = ±6ℓ and ξ 3 = 0, -6ℓ. This testing array, henceforth referred to as T 6 , is used for both source and sensor locations, so that Q = M = 602 in this case. For clarity, it is worth noting that the three-dimensional region in Fig. 2 is that enclosed by T [START_REF] Bonnet | Boundary integral equations methods for solids and fluids[END_REF] .

In what follows, the topological sensitivity is evaluated at the nodes of a uniform sampling grid G consisting of 51 × 51 = 2, 601 points in the 2D (sectional) diagrams and 51 × 51 × 51 = 132, 651 points in the 3D plots. The grid samples the region -5ℓ ≤ ξ 1 , ξ 2 ≤ 5ℓ, -5.5ℓ ≤ ξ 3 ≤ -0.5ℓ, i.e. most of the cubical enclosed by testing configuration T 6 . For the highest wavenumber considered (kℓ = 8), the horizontal and vertical spacing of the sampling grid are approximately 1/4 and 1/8 of the wavelength, respectively. Except for the results shown in Fig. 7, the topological sensitivity field is computed on the basis of a spherical nucleating inclusion shape B. All computations, including the BEM used for the creation of synthetic data, are implemented within the matlab environment and run on a laptop computer with a 1.7MHz Pentium IV processor. On this platform, computation of T f over the entire 3D sampling grid G (132,651 points) for testing configuration T 6 in the acoustic half-space case takes about 30 minutes of CPU time, with the BEM calculation of synthetic data requiring approximately 15 additional minutes for the dual-scatterer configurations.

Influence of testing grid aperture

Testing configuration T 6 , defined earlier, is deemed to be of full aperture inasmuch as the source and sensor locations completely surround the region of interest. To study the effects of limited aperture, two other testing configurations, T 1 and T 2 , are considered. Configuration T 1 consists of a uniform (10×10) square "source" grid of size 10ℓ with 11×11 source points, and a uniform (20×20) square "receiver" grid of size 10.5ℓ with 22×22 sensor locations, both grids lying in the plane ξ 3 = 0 and centered at

ξ 1 = ξ 2 = 0. Configuration T 2 = T 1 ∪T ′ 1 ,
where T ′ 1 is obtained by translating T 1 in the amount of 6ℓ in the negative ξ 3 -direction. Accordingly, (Q, M ) = (121, 484) for T 1 and (242, 968) for T 2 .

Figures 3 and4 shows the distribution of T f (x o , β true , γ true ) in the vertical plane (2, 0.5) (0, 1) (0, 1) Centroid (-2, 0.5 , -3) (2, 0.5 , -4) (-2, 0.5 , -3) (-2, 0.5 , -3) (-2, 0.5 , -3) x true c /ℓ (2, 0.5 , -4) (2, 0.5 , -4) (2, 0.5 , -4) Semi-axes (0.8, 0.4 , 0.4) (0.4, 0.4 , 0.8) (0.8, 0.4 , 0.4) (0.8, 0.4 , 0.4) (0.8, 0.4 , 0.4) a true /ℓ (0.4, 0.4 , 0.8) (0.4, 0.4 , 0.8) (0.4, 0.4 , 0.8) ξ 2 = 0.5ℓ, computed over the 51×51 uniform sampling grid for the testing arrangements T 1 , T 2 and T 6 , wavenumbers kℓ = 2, 4, 8, and scatterer configurations P-F (Fig. 3) and P-H (Fig. 4). It is not surprising that the best delineation of the obstacle for both configurations is obtained for T 6 , i.e. the case of full aperture. The differences between the images obtained under the respective conditions of full-space and halfspace background domains are highlighted by the detail in Fig. 5. In the display, the distribution of T f (x o , β true , γ true ) is seen to provide a good delineation of the scatterer shape for the full-space case, especially for kℓ = 8. In contrast, the distribution of T f (x o , β true , γ true ) for the half-space case is skewed towards the "top" surface ξ 3 = 0. On denoting the largest dimension of the scatterers by L = 1.6ℓ, one may note that the approximate wavelengths of "illuminating" waves featured in Figs. 3 and4 are 4L, 2L, L and 0.5L for kℓ = 1, 2, 4 and 8, respectively. In particular, the wavelengths for kℓ = 1, 2 and 4 are larger than or of a comparable size to the diameter of the scatterer and can thus be categorized as belonging to the so-called resonance region [START_REF] Colton | Inverse acoustic and electromagnetic scattering teory[END_REF].

Influence of testing grid density

With reference to the half-space scatterer configuration P-H and the full-aperture testing arrangement T 6 , the distributions of T f (x o , β true , γ true ) are computed for three densities of source and measurement grids, namely T 6 (2), T 6 (5) and T 6 (10) ≡ T 6 and plotted in Fig. 6. Here T 6 (2) and T 6 [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] are defined in the same way as T 6 where Q = M = 602 but on the basis of 2×2 (Q = M = 26) and 5×5 (Q = M = 152) rectangular grids on each face, respectively. From the display, one may note that the coarser grids T 6 (2) and T 6 (5) still provide acceptable results (with the notable exception of T 6 (2) for kℓ = 8) and in particular better than those obtained using "dense" but partial-aperture grids T 1 and T 2 , despite the larger number of source and sensor points in the latter two cases. These results suggest that the testing aperture is a more critical factor than the grid density for quality of obstacle delineation on the basis of topological sensitivity distribution.

Nucleating obstacle: influence of shape

In assessing the effectiveness of generalized topological sensitivity as a tool for preliminary obstacle reconstruction, a practical question arises as to the effect of the (assumed) shape B of a nucleating obstacle, manifested via tensor A in [START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF], on the spatial distribution of T f . This influence is illustrated in Fig. 7 which shows the distribution of T f (x o , β true , γ true ) in the vertical plane ξ 2 = 0.5ℓ, on the basis of the testing arrangement T 6 , for three reference shapes B: sphere (a 1 = a 2 = a 3 ), "vertical" ellipsoid (2a 1 = 2a 2 = a 3 ), and "horizontal" ellipsoid (a 1 = 2a 2 = 2a 3 ), where a 1 , a 2 and a 3 denote again the ellipsoid semi-axes. As can be seen from the display, the shape of B has a significant effect on the overall magnitude of T f , but not on its relative distribution. This observation is supported by the similar results obtained for other obstacle combinations, omitted here for brevity. Since it is the relative distribution of T f that provides most information about the support of the obstacle, one of the conclusions of this study is that the choice of B is not a critical component of the proposed (qualitative) imaging technique. This conclusion, however, may not hold for procedures where higher-order terms are considered in the expansion [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF], which may furnish more accurate information about the geometry of the scatterer (see [START_REF] Bonnet | Small-obstacle expansion in 3-d inverse scattering[END_REF] for impenetrable obstacles).

Nucleating obstacle: influence of constitutive parameters

In Figs. 34567, the T f -distribution is calculated assuming "correct" material parameters for the nucleating obstacle, namely (β, γ) = (β true , γ true ). The aim of this section is to examine the effect of material constants β and γ on the generalized topological sensitivity, and whether considering the T f -distributions with spatially-varying (β, γ) allows one to obtain information about the constitutive nature of the true scatterer.

With reference to the half-space scatterer configuration P-H, the distribution of T f (x o , β true , γ true ) is first compared in Fig. 8 to T f (x o , 0, γ true ) obtained on the basis of a rigid infinitesimal obstacle. As mentioned earlier, all simulations are performed on the basis of the testing arrangement T 6 = T 6 [START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF]. The comparison indicates that T f (x o , β true , γ true ) provides a markedly better indication of the obstacle location and shape, especially for kℓ = 4, 8.

The influence of β on T f is examined further in Fig. 9, where T f (x true c , β, γ true ) (evaluated at the centroid of the true obstacle) is plotted against β for scatterer configurations P-F, P-H, R-F, R-H and resonance-region wavenumbers kℓ = 1, 2, 4. Without exception, the band-limited "optimal" values β = β opt where T f (x true c ; β, γ true ) takes the maximum negative value are consistent with the mass densities of respective true scatterers (P: β true = 2, R:

β true = 0).
For completeness, the effect of γ on T f is examined in terms of T f (x true c , β true , γ), again evaluated at the centroid of the true scatterer. In this case, all diagrams are by default linear in terms of γ 2 , see [START_REF] Hobson | The theory of spherical and ellipsoidal harmonics[END_REF] or [START_REF] Volkov | Numerical methods for locating small dielectric inhomogeneities[END_REF]. Unfortunately, the synthetic results generated for a number of configurations indicate that the band-limited "optimal" value of γ is not always consistent with the constitutive nature of the respective true obstacle.

The foregoing results, and in particular those in Fig. 9, lend themselves to the idea of computing the "optimal" T f -distribution, obtained via point-wise minimization of T f (x o , β, γ) with respect to β so that

T f,opt (x o ) ≡ T f (x o , β opt (x o ), γ true ),
with the search interval for β opt taken here as 0 ≤ β ≤ 4. In this context, the key question is whether the "optimized" distribution of β opt (x o ) thus obtained is consistent with the true scatterer configuration. To facilitate the graphical interpretation, a thresholded variant of

T f,opt (x o ) is introduced according to Tf,opt (x o ) = T f,opt (x o ) (T f,opt ≤ CT min f,opt ), Tf,opt (x o ) = 0 (T f,opt > CT min f,opt ) (38)
where T min f,opt = min x o ∈G T f,opt (x o ) and C = 0.2 + (kℓ) -1 . In this way, an approximate image of the obstacle is formed through an assembly of sampling points where the optimal topological sensitivity T min f,opt takes "sufficiently" high negative values relative to the global minimum, T min f,opt < 0. The featured choice of a frequency-dependent threshold C is established on an ad-hoc basis, and reflects the fact that the T f -distributions are found to be more "smeared" at lower frequencies (see Figs. 3456). On the basis of (38), a corresponding (thresholded) distribution of the optimal β-parameter is taken as

βopt (x o ) = 0 ( Tf,opt (x o ) < 0), βopt (x o ) = 1 ( Tf,opt (x o ) = 0). ( 39
)
The distributions of Tf,opt and βopt , obtained for the single-scatterer configurations P-F and R-F, are shown in Figs. 10 and 11 for kℓ = 2, 4, 8. For the resonance-region wavenumbers (kℓ = 2, 4), the distributions of Tf,opt reasonably approximate the support of respective scatterers. The corresponding βopt -distributions likewise take values that are consistent with the qualitative nature of the respective "true" obstacles. The results for the "intermediate" wavenumber kℓ = 8, on the other hand, are mixed and are not deemed as effective as those obtained for the resonance region. The distributions of Tf,opt and βopt obtained for the dual-scatterer configurations PP-H, PR-H and RR-H, shown in Figs. 12-14, lead to the same conclusion. In particular, the kℓ = 2, 4 distributions of βopt obtained for configuration PR-H, which features one soft and one rigid obstacle, point to the correct nature of each scatterer.

For completeness, the level surfaces of T f,opt (x o ) defined by T f,opt = CT min f,opt where C = 0.2 + (kℓ) -1 are computed over the entire 3D sampling grid G and graphically shown in Fig. 15 for kℓ = 2, 8 and scatterer configurations PP-H, PR-H and RR-H. These three-dimensional plots in particular demonstrate that the values of topological sensitivity deemed significant (i.e. lower than the threshold value CT min f,opt ) do not occur outside a neighborhood of the scatterers, an observation which is found to hold for all configurations examined in this study.

Conclusions

In this paper, the concept of topological sensitivity is generalized to permit nucleation of penetrable obstacles as it pertains to the three-dimensional Helmholtz equation and inverse scattering problems in acoustics. In the context of arbitrarily-shaped inclusions, the featured formula is shown to consist of a dipole and a monopole term, related respectively to the density and compressibility contrast between the nucleating obstacle and a background medium. Explicit formulas are given for special cases when the nucleating inclusion takes spherical or ellipsoidal shape. The introduction of an adjoint solution further permits to consider nucleation of arbitrarily-shaped inclusions in an infinite, semi-infinite or finite background medium. To highlight the utility of proposed developments in dealing with inverse scattering problems, a set of numerical results is included wherein hidden obstacles are exposed through regions where (closedform) topological sensitivity attains negative values. On varying the mass density of a nucleating obstacle, it is also shown that the proposed methodology can be used in conjunction with long wavelengths (the so-called resonance region) for both geometric and material identification.

T 1 T 2 T 6 kℓ = 2 -5 -4 -3 -2 -1 0 1 2 3
x 10 (ℓρ

-3 -5 -4 -3 -2 -1 0 1 2 3 4 5 -10.5 -10 -9.5 -9 -8.5 -8 -7.5 -7 -6.5 -6 -5.5
-2 c -4 ) T f (x o , β true , γ true ) (ℓρ -2 c -4 ) T f (x o , 0, γ true ) kℓ = 2 -5 -4 -3 -2 -1 0 1 2 3

Figure 1 .

 1 Figure 1. Configuration of the scattering problem

Figure 2 .

 2 Figure 2. True obstacle configurations P (a), R (b), and PR (a-b).

Figure 3 .

 3 Figure 3. Influence of testing aperture: distribution of (ℓρ -2 c -4 ) T f (x o , β true , γ true ) in the vertical plane ξ 2 = 0.5ℓ for the scatterer configuration P-F.

Figure 4 .

 4 Figure 4. Influence of testing aperture: distribution of (ℓρ -2 c -4 ) T f (x o , β true , γ true ) in the vertical plane ξ 2 = 0.5ℓ for the scatterer configuration P-H.

T 6

 6

Figure 6 .Figure 7 .

 67 Figure 6. Influence of testing grid density (P-H configuration): distribution of (ℓρ -2 c -4 ) T f (x o , β true , γ true ) in the vertical plane ξ 2 = 0.5ℓ, calculated for the testing arrangements T 6 (2), T 6 (5) and T 6 (10) = T 6 .

Figure 8 .Figure 9 .

 89 Figure 8. Influence of trial obstacle parameters (P-H configuration): distribution of T f (x o , β true , γ true ) (left) and T f (x o , 0, γ true ) (right) in the vertical plane ξ 2 = 0.5ℓ.

(ℓρ - 2 cFigure 10 .

 210 Figure 10. Distribution of Tf,opt (x o ) and β opt (x o ) (P-F configuration).

(ℓρ - 2 c - 4 )Figure 11 .

 2411 Figure 11. Distribution of Tf,opt (x o ) and β opt (x o ) (R-F configuration).

(ℓρ - 2 c - 4 )Figure 12 .

 2412 Figure 12. Distribution of Tf,opt (x o ) and β opt (x o ) (PP-H configuration).

(ℓρ - 2 c - 4 )Figure 13 .

 2413 Figure 13. Distribution of Tf,opt (x o ) and β opt (x o ) (PR-H configuration).

(ℓρ - 2 c - 4 )Figure 14 .

 2414 Figure 14. Distribution of Tf,opt (x o ) and β opt (x o ) (RR-H configuration).

Figure 15 .

 15 Figure 15. Level surfaces of T f,opt (x o ) defined by T f,opt = C(kℓ)T min f,opt , for kℓ = 2, 8 and scatterer configurations PP-H, PR-H and RR-H.

Table 1 .

 1 Scatterer configurations.

	Name	P	R	PP	RR	PR
	Properties	(2, 0.5)	(0, 1)	(2, 0.5)	(0, 1)	(2, 0.5)
	(β true , γ true )

† Inverse Problems 22:1761-1785 (2006)

Acknowledgments

The support provided by the National Science Foundation through grants CMS-0324348 and INT-0340590 to B. Guzina and the National Center for Scientific Research, France's Ministry of Research, through award No. 17151 to M. Bonnet is gratefully acknowledged. The authors are also thankful to the University of Minnesota Supercomputing Institute for the computing resources made available during the course of this investigation.


Sylvain Takerkart

Guillaume Auzias

Lucile Brun

Olivier Coulon

Mapping Cortical Shape Differences Using a Searchlight Approach Based On Classification of Sulcal Pit Graphs

Keywords: brain, mri, classification, mapping, searchlight, cortex, anatomy, sulcus, sulcal pit, graph, kernel

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In the past few years, the topography of the cortical surface has raised a lot of interest, in particular to find biomarkers of pathologies [START_REF] Im | Quantification and Discrimination of Abnormal Sulcal Patterns in Polymicrogyria[END_REF][START_REF] Auzias | Atypical sulcal anatomy in young children with autism spectrum disorder[END_REF] or exhibit differences associated to functional specificities [START_REF] Sun | The effect of handedness on the shape of the central sulcus[END_REF]. Behind the large apparent variability of cortical folding patterns, a specific attention has been brought to the deepest part of sulci, either to elaborate theoretical models of cortical anatomy and development [START_REF] Régis | Sulcal root generic model: a hypothesis to overcome the variability of the human cortex folding patterns[END_REF], or to automatically extract robust cortical landmarks. For the latter, the work of Im has been particularly important in defining the concept of sulcal pits. Sulcal pits are defined as local maxima of depth within each cortical fold [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF]. They can be extracted, together with their associated sulcal basins, via a watershed algorithm performed on a sulcal depth map defined on the mesh of the cortical surface [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF][START_REF] Auzias | Deep sulcal landmarks: algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF]. Sulcal pits have ben linked with genetic factors [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF] or developmental pathologies [START_REF] Im | Quantification and Discrimination of Abnormal Sulcal Patterns in Polymicrogyria[END_REF].

The challenge in using sulcal pits is to find stable patterns within a population, despite the apparent inter-subject variability. To do so, it is essential to have a good representation of folding patterns. Sulcal pits define sulcal basins that parcellate the surface of the neo-cortex (Fig. 1), and the adjacency of neighboring sulcal basins can be used to define sulcal pit graphs. Such representations of local anatomical patterns have been successfully analysed at the brain lobe level using a spectral graph-matching technique [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF][START_REF] Im | Quantification and Discrimination of Abnormal Sulcal Patterns in Polymicrogyria[END_REF].

In this work we introduce a multivariate technique that enables to localize more precisely the differences between populations from local patterns of sulcal pits. It relies on two main contributions: i) the design of graph kernel that allows the group classification to be performed directly in graph space; this kernel has very few parameters that can be efficientely inferred from the data; ii) the definition of a searchlight mapping strategy [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF] that enables to localize the discriminant patterns. In what follows, we describe our method and demonstrate its effectiveness by mapping differences between male and female subjects using a large dataset.

MATERIALS AND METHODS

Extracting sulcal pits

In order to localize the sulcal pits, we used a modified version of the procedure initially proposed in [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF], designed to yield reproducible sulcal pits [START_REF] Auzias | Deep sulcal landmarks: algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF]. First, we estimated the sulcal depth map for each subject using the depth potential function, a measure that integrates both curvature and convexity information [START_REF] Boucher | Depth potential function for folding pattern representation, registration and analysis[END_REF]. Then, we applied a watershed by flooding algorithm to extract depth maxima and their corresponding sulcal basins on the mesh (Fig. 1). A merging of sulcal basins was performed during the flooding in order to filter spurious extrema caused by noise in the depth map. Merging parameters were optimized as described in [START_REF] Auzias | Deep sulcal landmarks: algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF] (see also [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF]).

For each subject, we obtain a set of pits P = {p i } i=N i=1 and their basins {B i } i=N i=1 (∀i, p i is the deepest point of B i ). Note that i) the basins form a complete parcellation of the cortical surface and ii) the number of pits N can vary across subjects.

Representing patterns of sulcal pits as graphs

A natural way to formally represent a pattern formed by a set of pits, including their spatial organization, is to construct a graph. We here use the method proposed in [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF]. For a given subset P ⊂ P of M pits (M ≤ N ), we define a node of the graph for each pit p i ∈ P. The graph edges are then given by the spatial adjacency of their associated basins: this defines a binary adjacency matrix A = (a ij ) ∈ R M ×M (a ij = 1 if B i and B j are adjacent, and 0 otherwise), that encodes the spatial organization of neighboring pits. In order to better characterize the pattern of pits, we add two attributes to each graph node: i) its depth d i , because it is an intrinsic characteristic of the pit, and ii) the 3d coordinates X i of the corresponding vertex on the sphere after surface basedalignment to a common template performed in freesurfer, so that we can compare the locations of the graph nodes across subjects. Let D = {d i } ∈ R M be the vector of depth values and X = {X i } ∈ R M ×3 be the matrix of coordinates of all graph nodes.

Any pattern of sulcal pits can therefore be fully represented by an attributed graph defined as G = ( P, A, D, X ). Examples of sulcal pit graphs are shown on Fig. 1.

Graph-based support vector classification

Our first contribution is to introduce a new graph kernel that measures the similarity between two graphs

G = ( PG , A G , D G , X G ) and H = ( PH , A H , D H , X H ).
Equipped with such a kernel, one can use the kernel trick to perform classification directly in graph space using a support vector machine (see for instance [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]). Similarly to previous work dedicated to inter-subject multi-voxel pattern analysis of functional MRI [START_REF] Takerkart | Graph-based inter-subject pattern analysis of fMRI data[END_REF], we designed a kernel that exploits all the features of the sulcal pit graphs defined in 2.2.

Kernel definition. Our kernel compares all pairs of nodes (i.e potential edges) g ij and h kl , respectively in G and H, without trying to perform a one-to-one pit-matching as in [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF]. As such, it belongs to the class of walk-based graph kernels [START_REF] Gärtner | On graph kernels: Hardness results and efficient alternatives[END_REF] and uses the most elementary walks, of length one. It combines (see eq. 4) the different features of the graphs by using several sub-kernels within the convolution kernel framework [START_REF] Haussler | Convolution kernels on discrete structures[END_REF]. A first sub-kernel aims at ensuring that the comparisons are performed only if g ij and h kl are actual edges. This is done with the linear kernel on the binary entries of the adjacency matrices:

k a (g ij , h kl ) = a G ij .a H kl , (1)
which takes the value 1 if a G ij = a H kl = 1 and 0 otherwise. A second sub-kernel uses a product of gaussian kernels on the coordinates of the nodes of g ij and h kl :

k x (g ij , h kl ) = e -X G i -X H k 2 /2σ 2 x • e -X G j -X H l 2 /2σ 2 x . (2)
In practice, k x acts as a spatial filter that weights the comparisons of edges with their proximity, thus eliminating the comparisons of edges that are far away from eachother and allowing for inter-subject variability (if edges are close, but not perfectly matched across subjects). Finally, the last subkernel compares the depth attributes using the same principle:

k d (g ij , h kl ) = e -d G i -d H k 2 /2σ 2 d • e -d G j -d H l 2 /2σ 2 d (3)
The full kernel is defined as the combination of the three sub-kernels applied on all pairs of nodes of G and H:

K(G, H) = M G i,j=1 M H k,l=1 k a (g ij , h kl )•k x (g ij , h kl )•k d (g ij , h kl ) (4)
Note that the number of nodes M G and M H in G and H can be different.

Estimating hyper-parameters. The value of the two bandwidth parameters σ x and σ d need to be chosen. We adapt a standard heuristic (used for vector inputs to a gaussian kernel) that selects the median euclidean distance between all observations in the training dataset. We here use for σ x and σ d the median euclidean distance between the coordinates and depth attributes of the nodes in all graphs available at training time. The values of these parameters are therefore estimated directly and easily from the training dataset, using this heuristic which was proved efficient [START_REF] Takerkart | Graph-based inter-subject pattern analysis of fMRI data[END_REF].

Searchlight mapping

Our second contribution consists in setting up a searchlight method aimed at mapping the discriminative power of local patterns formed by sulcal pits across populations. The socalled searchlight method [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF] consists in using a multivariate statistical model in a sliding window that defines a local neighborhood and moves along the space of interest. A summary statistic computed from this model (for instance the accuracy of a classifier) is then assigned to the center of the window, thus yielding a spatial map that allows the localization of the informative patterns. More specifically, in order to construct a searchlight scheme for a given task, one needs to define i) the space of interest, ii) the spatial sampling strategy of this space, iii) how to define local patterns at each location, iv) the statistical model that addresses the task itself (in our case: group discrimination) and v) the summary statistic to be mapped onto the space of interest. Below, we instantiate these five items to define our pit-based searchlight method.

Since sulcal pits are defined as local depth extrema on the cortex, it is natural to perform our searchlight strategy along the cortical surface (i). Because the size of the sulcal patterns of interest are far larger than the triangles comprising the cortical mesh, it is not necessary to sample the cortex using each vertex of the mesh. We therefore consider a coarser set of Q evenly distributed points over the spherical template provided with freesurfer and map those onto the cortex of each subject (ii). At each of these points and for each subject, we apply the graph construction described in 2.2 using the set of pits P located within a given radius R from the center of the neigborhood (iii). Then, we estimate a non linear support vector classifier to discrimate the populations based on the local sulcal pit graphs, using our graph kernel 2.3 (iv). We compute the classifier accuracy using a cross-validation, and use the average accuracy ā across folds as the summary statistic to build our information map (v). In order to assess whether this accuracy is significantely above chance level, we use a non parametric permutation test that provides a p-value for the null hypothesis of no differences between populations.

EXPERIMENTAL RESULTS AND DISCUSSION

Mapping gender differences

The Open Access Series of Imaging Studies (OASIS, www. oasis-brains.org) cross-sectional database offers a collection of 416 subjects aged from 18 to 96. For each subject, three to four individual T1-weighted MP-RAGE scans were obtained on a 1.5T Vision system (Siemens, Erlangen, Germany) with the following protocol: in-plane resolution = 256x256 (1 mm x 1 mm), slice thickness = 1.25 mm, TR = 9.7 ms, TE = 4 ms, flip angle = 10, TI = 20 ms, TD = 200 ms. Images were co-registered and averaged to create a single image with a high contrast-to-noise ratio. From this database, we selected two groups of 67 male and 67 female healthy right-handed subjects, aged 18 to 34, matched in age.

Using this data, our goal was to map gender differences from local patterns of sulcal pits. At each point of the searchlight, a binary classification problem (males vs. females) was solved using the graph-based classifier defined previously. The average classification accuracy ā was estimated using a 10-fold cross-validation that ensured gender balance in both the training and test datasets for each fold. The number of searchlight locations Q on the cortex was set to 300, and the radius of the neighborhood to R = 50mm.

Results and discussion

We found several regions for which the classification accuracy ā was significantely above chance level (0.5), i.e for which the probability of the null hypothesis of no difference between males and females was p < 0.001. The corresponding accuracy scores were between ā = 0.61 and ā = 0.63. As visible on Fig. 2, the three main regions are : i) a compact set of 4 searchlight locations at the dorsal anterior extremity of the right frontal lobe, ii) a location at the junction between the right inferior temporal sulcus and the occipital-temporal junction, and iii) a location in the left calcarine sulcus. Although previous studies using standard voxel-or surface-based morphometry [START_REF] Ruigrok | A meta-analysis of sex differences in human brain structure[END_REF][START_REF] Escorial | Sex differences in neocortical structure and cognitive performance: A Surface-Based Morphometry Study[END_REF] do not report differences directly related to folding, they indicate that several other anatomical characteristics (grey matter volume, density, thickness or cortical area) correlate with our results by finding corresponding regions. In particular, local differences in grey matter volume are to some extent predictive of the differences shown by our results, because they embed information about both cortical surface area, which is linked to folding, and cortical thickness, which can induce local depth changes. Therefore, our searchlight scheme is able to detect local anatomical differences between genders that are consistent with the standard morphometry literature. However, another study that examined group differences between sulcal pit graphs did not find any differences between genders (auxiliary result in [START_REF] Im | Quantification and Discrimination of Abnormal Sulcal Patterns in Polymicrogyria[END_REF]). We list hereafter four potential reasons that could explain this disagreement: i) our dataset was much larger (134 vs. 26 subjects); ii) our pit detection algorithm is slightly different [START_REF] Auzias | Deep sulcal landmarks: algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF]; iii) our kernel provides a graph similarity measure that might be sensitive to distinct properties of the graphs; iv) their study focused on brain lobes, whereas we consider smaller neighborhoods; it might be that the discriminant information lies at a finer scale than brain lobes. Further investigations will be needed to decipher which of these four factors contribute to the gain in sensitivity apparently offered by our framework. This is the first time that a searchlight approach is proposed for studying local patterns of cortical shape. The main strength of this approach is to overcome the intrinsic limitations of ROI-based approaches since the regular sampling of the surface implies an overlap between patterns. Furthermore, varying the size of the neighborhood (here the radius R) enables to study different spatial scales, which is necessary in clinical applications to deal with various pathologies. The second strength of our approach is the use of a similarity measure that is a positive definite kernel. This opens to using any kernel method, meaning we can tackle various tasks such as classification but also regression or clustering. Our approach is therefore powerful for localizing effects and versatile because opened to numerous applications.

CONCLUSION

We have introduced a multivariate searchlight mapping technique that uses a graph-based support vector classifier to localize informative patterns of sulcal pits. We have showed that our framework can detect anatomical gender differences that are consistent with the literature. This versatile and powerful pit-based searchlight approach should therefore find numerous applications related to clinical populations.

Fig. 1 .

 1 Fig. 1. A. Illustration of the sulcal pits and their basins on one subject. B. Examples of sulcal pit graphs in the right frontal lobe for four subjects (the color of the nodes encodes the relative depth of the pit, from red (deep) to blue (shallow).

Fig. 2 .

 2 Fig. 2. Mapping of the informative patterns of sulcal pits for gender differences. Each sphere is one of the Q = 300 searchlight locations. Their color encodes the classification accuracy ā. The large spheres indicate the locations for which ā is significantely above chance level (p < 0.001).