## Abstract

This study deciphers the topological sensitivity (TS) as a tool for the reconstruction and characterization of impenetrable anomalies in the high-frequency regime. It is assumed that the anomaly is simply connected and convex, and that the measurements of the scattered field are of the far-field type. In this setting, the formula for TS—which quantifies the perturbation of a cost functional due to a point-like impenetrable scatterer—is expressed as a pair of nested surface integrals: one taken over the boundary of a hidden obstacle, and the other over the measurement surface. Using multipole expansion, the latter integral is reduced to a set of antilinear forms featuring Green's function and its gradient. The remaining expression is distilled by evaluating the scattered field on the surface of an obstacle via Kirchhoff approximation, and pursuing an asymptotic expansion of the resulting Fourier integral. In this way, the TS is found to survive upon three asymptotic lynchpins, namely (i) the near-boundary approximation for sampling points close to the ‘exposed’ surface of an obstacle; (ii) uniform expansions synthesizing the diffraction catastrophes for sampling points near caustic surfaces, lines and points; and (iii) stationary phase approximation. Within the framework of catastrophe theory, it is shown that, in the case of the full source aperture, the TS is asymptotically dominated by the (explicit) near-boundary term—which explains the previously reported reconstruction capabilities of this class of indicator functionals. The analysis further shows that, when the (Dirichlet or Neumann) character of an anomaly is unknown beforehand, the latter can be effectively exposed by assuming point-like Dirichlet perturbation and considering the sign of the leading term inside the reconstruction.

## 1. Introduction

Waveform tomography and in particular inverse obstacle scattering are essential to a broad spectrum of science and technology disciplines, including geophysics, oceanography, optics, aeronautics and non-destructive material testing. In general, the relationship between the wavefield scattered by an obstacle and its geometry (or physical characteristics) is nonlinear, which invites two overt solution strategies: (i) linearization via, e.g. Born approximation and ray theory [1] or (ii) pursuit of the nonlinear minimization approach [2]. Over the past two decades, however, a number of *sampling methods* have emerged that both consider the nonlinear nature of the inverse problem and dispense with iterations. Commonly, these techniques deploy an indicator functional that varies with spatial coordinates of the trial, i.e. sampling point, and projects observations of the scattered field onto a functional space reflecting the ‘baseline’ wave motion in a background medium. This indicator functional, designed to reach extreme values when the sampling point strikes the anomaly, accordingly provides a tomogram via its (thresholded) spatial distribution. Examples of such imaging paradigm include the linear sampling method [3] and the factorization technique [4] in the context of extended (i.e. finite-sized) scatterers, as well as the MUSIC algorithm [5] and the direct approach [6] as techniques catering primarily for point-like targets.

Another sampling take on inverse scattering, which motivates this study, is the method of topological sensitivity (TS) [7,8]. Stemming from the framework of shape optimization [9], this technique has emerged as an effective tool for the waveform tomography of extended obstacles in acoustics [8,10–13], electromagnetism [14,15] and elastodynamics [16–19]. Formally, the TS quantifies the leading-order perturbation of a given misfit functional when an infinitesimal scatterer is introduced at a sampling point in the reference domain—being imaged for obstacles. From the application viewpoint, the appeal of TS resides in its forthright computability as a bilinear form in terms of two (free and adjoint) forward solutions for the reference domain. Following the heuristic argument, this quantity is then used as *obstacle indicator* by identifying the support of its pronounced negative values with an anomaly.

Despite the mounting numerical [8,15–18] and experimental [20,21] evidence of the imaging capability of the method in a variety of sensing configurations, a theoretical justification of the TS as an obstacle indicator function is still lacking. So far, Dominguez *et al*. [22] established the analogy between the TS and time reversal, while Bellis *et al*. [23] elucidated the TS reconstruction of (i) point-like anomalies and (ii) extended weak anomalies in the sense of small material contrast (Born approximation) and/or low excitation frequency. Further, Ammari *et al.* [24] explained how the TS discerns small acoustic obstacles, and in [25] they exposed the link between the TS and error backpropagation. To date, however, the reported ability of TS to reconstruct obstacles of arbitrary (finite) size and contrast has eluded both physical understanding and rigorous justification. The problem is highlighted by the repeated observations [8,11,17,21] that at higher frequencies, the usual reconstruction heuristics does not apply for the negative TS values tend to localize in a narrow region ‘about the boundary’ of an anomaly [11]—rather than canvassing its support.

To help bridge the gap, this study focuses on the imaging by TS of impenetrable (Dirichlet or Neumann) obstacles in the short-wavelength regime. First, the expression for the germane indicator functional is reduced, via multipole expansion and Kirchhoff approximation, to a Fourier-type surface integral over the illuminated part of the anomaly's boundary. Making use of the high-wavenumber hypothesis, the latter is distilled to three distinct asymptotic representations, namely (i) the near-boundary approximation for sampling points within few wavelengths from the ‘exposed’ surface of an obstacle; (ii) diffraction catastrophes (of codimension less than four) for sampling points near caustic surfaces, lines and points; and (iii) non-uniform, i.e. stationary phase approximation. Under the premise of a single illuminating (plane) wave, it is found that the distribution of TS, while carrying hints about the shape of an anomaly via the near-boundary contribution, is controlled by the caustics. By the way of the catastrophe theory and diffraction scaling laws, on the other hand, it is shown that in the case of the *full aperture* of illuminating wavefields, the TS is asymptotically dominated by the (explicit) near-boundary term—which explains the evidenced imaging capabilities of this class of indicator functionals. This result further unveils the new reconstruction logic at short wavelengths, where the *boundary* of an anomaly is obtained as a *zero level set* of the TS field separating its extreme negative and extreme positive values. From the practical point of view, such paradigm allows for obstacle reconstruction without the use of an ad-hoc threshold parameter, where the sign of TS inside the reconstruction can be used to identify the anomaly type (if unknown) as either Neumann or Dirichlet. The analysis is accompanied by numerical results and an application towards obstacle reconstruction to a recent set of experimental data [21].

## 2. Preliminaries

Consider the inverse scattering of time-harmonic scalar waves by a simply connected, convex and impenetrable obstacle *S*=∂*D*, where *R*_{1} centred at the origin. On denoting by *u*^{i} on *D*, the total field
*R*_{2}=*α*^{−1}*R*_{1} (*α*<1) centred at the origin (figure 1). The reference background medium is assumed to be homogeneous with wave speed *c* and mass density *ρ*.

*Objective.* The goal is to reconstruct *D*, and to identify its character (as either sound-soft or sound-hard) in situations when such prior information is unavailable.

*Sensory data.* Writing the implicit time dependence as e^{iωt}, the incident field is assumed in the form of a *plane wave*, *u*^{i}=e^{−ikξ⋅d}, endowed with wavenumber *k*=*ω*/*c* and direction ** d**∈

*Ω*, where

*Ω*is a unit sphere. For each

**, values of the total field**

*d**u*are collected over

*Γ*

^{obs}.

*Cost functional.* To help solve the inverse problem, consider the misfit functional
** d**, where

*trial*(Dirichlet or Neumann) obstacle;

*v*is the total field generated by the action of

*u*

^{i}on

*φ*is a distance function that is differentiable with respect to the real and imaginary parts of its first argument. In what follows,

*φ*is assumed to take the usual least-squares format

*Green's function.* For further reference, let
*k*, where ** x** signifies the source location;

*r*=|

**−**

*ξ***|, and ∇**

*x**G*indicates differentiation with respect to

*the first*argument.

*Dimensional platform.* In the sequel, all quantities are assumed to be dimensionless. This is accomplished by taking the radius of the inner sphere, the mass density of the background medium and the sound speed in the background medium as the reference length, mass density and velocity. In this setting, one in particular has *R*_{1}=1 and *R*_{2}=*α*^{−1}.

### (a) Topological sensitivity

Let *J*(∅) due to insertion of a vanishing impenetrable obstacle *T*(** x**°), the so-called TS, is independent of

*ϵ*. Taking

*unit ball*, analyses further show [6,10,26,27] that

*A*and

*B*depend on the physical character of a vanishing perturbation, and

*u*

^{a}is the so-called adjoint field, generated in

*Γ*

^{obs}. Note that the implicit assumption underpinning (2.5) is that the wavenumber

*k*is fixed in the limiting process, whereby

*kϵ*→0 in (2.4) [28,29]. As a prelude to the ensuing discussion, table 1 specifies the coefficients

*A*and

*B*in situations when

### (b) Scaled topological sensitivity indicator

As indicated earlier, the goal of this work is to provide an overarching high-frequency treatment of impenetrable (Dirichlet and Neumann) obstacles via the concept of TS. Given the dependence of *B* on *k* in the Neumann case (table 1), however, one may conveniently introduce a *scaled counterpart* of (2.5), namely
*γ* and normalized expansion coefficients A and B are given in table 1. Relative to (2.5), the use of (2.6) is advantageous in that (i) the new expansion coefficients are *k*-independent and (ii) the asymptotic order of T, in terms of powers of *k*, is invariant with respect to the (Dirichlet or Neumann) character of the vanishing perturbation. In the remainder of this work, (2.6) is *used as the basis* for the high-frequency reconstruction and identification of impenetrable obstacles. Note that the featured scaling of TS by *k*^{γ} (besides simplifying the analysis) does not affect the inverse solution, as the anomalies are usually identified from the relative variation of TS (e.g. [8,10–13]).

To help expose the high-frequency behaviour of T, one may expand the adjoint field in (2.6) as
** n** is the unit

*outward*normal on

*S*;

*u*

_{,n}=

**⋅∇**

*n**u*, and

*G*(

**,**

*x***)=**

*y**G*(

**,**

*y***) and ∇**

*x**G*(

**,**

*x***)=−∇**

*y**G*(

**,**

*y***).**

*x*### (c) Approximation of the component integrals over *Γ*^{obs}

The purpose of this section is to reduce the TS formula (2.8) to a single, Fourier-type surface integral that is amenable to short-wavelength approximation. Note that analogous derivations can be found in [24,30]. From (2.3), it follows that
_{1}) and (2.9_{2}), respectively, by ** a**,

**)=(**

*b***,**

*ζ***°) from its companion yields**

*x**R*

_{1}=1 and

*R*

_{2}=

*α*

^{−1}>1, it can be shown via triangle inequality that |

*E*|<

*α*

^{2}/(

*k*(1−

*α*

^{2}))+

*α*

^{2}/2+

*O*(

*α*

^{4}). When

*k*≥

*O*(1), (2.11) accordingly yields the Helmholtz–Kirchhoff identity

*O*(

*α*

^{n}) residual.

On differentiating (2.11) with respect to ** x**° and

**, it can be similarly shown that**

*ζ**G*indicates differentiation with respect to the first argument,

*r*=|

**°−**

*x***| and**

*ζ***is the second-order identity tensor.**

*I*

## 3. High-frequency behaviour of topological sensitivity

As examined earlier, the objective of this work is to understand the previously reported high-frequency patterns of the TS indicator functional (e.g. [8,11]), obtained when (2.7) is applied to the scattered field data at wavelengths (2*π*/*k*) that are smaller than the characteristic size of an obstacle, *L*_{o}. In this vein, the ensuing analysis assumes a *separation of scales* in that *ϵ*≪2*π*/*k*≪*L*_{o}, where *ϵ*→0 is the size of a vanishing perturbation in (2.4) while 2*π*/*k*, however small, is *fixed* in the limiting process. With such premise, consider the scattering of a plane wave, *u*^{i}=*M* e^{−ikx⋅d}, by convex impenetrable obstacle *D* as in §2. Next, let ** n** signify the outward normal on

*S*=∂

*D*; let

*S*

^{f}(

**)={**

*d***∈**

*x**S*:

**(**

*n***)⋅**

*x***<0} be the ‘front’ (i.e. illuminated) part of**

*d**S*, and denote by

*S*

^{b}(

**)={**

*d***∈**

*x**S*:

**(**

*n***)⋅**

*x***≥0} its ‘back’ side.**

*d*Here it is noted that the high-frequency analysis of integrals similar to those featured in §2a has been recently proposed in [31] towards developing an iterative scheme for the multi-static imaging of extended penetrable targets. There, it is shown that the high-frequency data can be used to construct a good initial guess for the illuminated part of the inclusion.

*Dirichlet obstacle as a testbed.* To provide specificity for the analysis, it is hereon assumed that the hidden, i.e. extended obstacle *D* is *sound-soft*. The case of a Neumann (sound-hard) obstacle is addressed separately in §4d and electronic supplementary material, appendix E. As it turns out, however, the latter developments draw heavily from the Dirichlet analysis—and in fact require only a minimal amount of additional deliberation.

### Remark 3.1

When *D* is of Dirichlet type, it is natural to let the trial vanishing obstacle *D* (as an impenetrable anomaly) constitutes prior information that may not be available. As a result, it is of interest to proceed with the inverse scattering of sound-soft *D* while allowing (A,B) to assume either pair of values given in table 1. As will be shown in §4e, such paradigm paves the way towards a simultaneous reconstruction and characterization of extended impenetrable obstacles in situations when their character is unknown beforehand.

### (a) Kirchhoff approximation

When the extended obstacle is sound-soft and *kL*_{o}≫1, the physical optics (Kirchhoff) approximation [32] states that
*u*^{i}=e^{−ikx⋅d} and substituting (2.12) and (2.13) into (3.2), one finds that
*O*(*α*^{2}) approximation error stemming from (2.12) and (2.13) is tacit.

In what follows, the high-frequency behaviour of (3.3) is characterized explicitly, assuming both (i) illumination by a single incident wave and (ii) full source aperture when ** d**∈

*Ω*. Specifically, the analysis shows that in the latter case, the TS distribution is approximated by a

*closed-form expression*(the main contribution of this work) whose extreme values are localized in a neighbourhood of ∂

*D*, as suggested by numerical investigations [8,11]. This result is stated in theorem 4.10 for sound-soft obstacles, and in theorem 4.12 for sound-hard obstacles. For completeness, electronic supplementary material, appendix E examines the ramifications of assuming Kirchhoff approximation (3.1) on the claim of theorem 4.10.

### (b) Contribution of non-degenerate stationary points

Consider first the high-frequency behaviour of (3.3) when the sampling point ** x**° straddles the region of interest

*excluding*a ‘thin-shell’ neighbourhood of

*S*

^{f}, namely

*ϵ*=

*O*(

*k*

^{−1}) is a length scale to be specified later. In this setting, the analysis can be facilitated by recalling (2.3) and rewriting (3.3) as

*S*

^{f}in terms of curvilinear surface coordinates (

*η*

^{1},

*η*

^{2}) as

*g*

_{pq}are the covariant components of the metric tensor.

As examined in [32], the leading-order asymptotic behaviour of (3.6) for large *k* is governed by the nature of the integrand in the neighbourhood of three types of *critical points*, namely: (i) the stationary points on *S*^{f} where ∇_{η}(** ζ**⋅

**±**

*d**r*) vanishes; (ii) the points on

*S*

^{f}where the integrand fails to be differentiable; and (iii) all points on the closed curve ∂

*S*

^{f}—the boundary of

*S*

^{f}. By way of (3.5),

*r*≥

*ϵ*>0 whereby the integrands in (3.6) are differentiable everywhere. One may also note that the latter

*vanish*on ∂

*S*

^{f}due to multiplier

**⋅**

*d***. Following the analysis in [32], the leading contribution of ∂**

*n**S*

^{f}to

*J*

_{1}and

*J*

_{2}can accordingly be shown to behave as

*O*(

*k*

^{−2}) when

*k*is large. By contrast, the contribution of a non-degenerate stationary point

***∈**

*ζ**S*

^{f}to a two-dimensional Fourier integral

*non-uniform*asymptotic approximation (e.g. [33]) as

*O*(

*k*

^{−1}), namely

*A*

_{pq}∈{−2,0,2} is the difference between the numbers of positive and negative eigenvalues of

*A*

_{pq}. Accordingly the portion of (3.6) due to non-degenerate stationary points can be computed, to the leading order, by summing the contributions of type (3.8).

#### (i) Stationary points of (3.7)

To evaluate (3.6) via the method of stationary phase [32], it is noted that
*ζ*^{±}∈*S*^{f} the stationary point of e^{ik(ζ⋅d±r)}, this implies that *S*^{f}. Making use of the inequality ** d**⋅

**<0, one finds from (3.10) that**

*n**J*

_{1}and

*J*

_{2}feature two types of stationary points, namely

*uniquely*determined by the projection of

**° along**

*x***on**

*d**S*

^{f}. In the light of the implicit specification of

*ζ*^{±}

_{II}, on the other hand, integrals

*J*

_{1}and

*J*

_{2}may have

*multiple*stationary points of type II. To provide further insight into (3.11), let

*given*boundary point

**∈**

*ζ**S*

^{f}is the stationary point of (3.6). This is illustrated in figure 2 which shows that the I

^{−}and II

^{+}loci emanate from

*S*

^{f}towards the exterior of

*D*, while their I

^{+}and II

^{−}counterparts extend (initially) from

*S*

^{f}towards the interior of

*D*. One also may note that at the ‘apex’ of

*S*

^{f}, where

**=−**

*n***, locus I**

*d*^{−}(resp. I

^{+}) coincides with locus II

^{+}(resp. II

^{−}). Such coalescence, however, does not pose special problems since each of the component integrals in (3.6) will have a stationary point of

*either*type I or type II that in this case coincides with the apex of

*S*

^{f}.

*Stationary point of type I.* Recalling (3.12), the asymptotic behaviours of *J*_{1} and *J*_{2} entail the contribution of a unique stationary point *ζ*^{±}_{I} when ^{±} otherwise. The results in the electronic supplementary material, appendix A(a) show that in the former case
*ζ*^{±}_{I} to T as
*r* is separated from zero thanks to (3.5).

*Stationary point of type II ^{+}.* From the analysis in the electronic supplementary material, appendices A and C, one finds that

*J*

_{1}and

*J*

_{2}feature a

*unique*stationary point

^{+}otherwise. In particular, it is shown that

*ρ*

_{1/2}are the principal radii of curvature of

*S*

^{f}at

*ζ*^{+}

_{II}, while the roots

*r*

_{1}≥

*r*

_{2}>0 solve (A5) (see the electronic supplementary material, appendix A(a)). On the basis of (3.3), (3.6), (3.8) and (3.16), the contribution of

*ζ*^{+}

_{II}to T reads

^{±}do not contribute to the leading asymptotic behaviour of TS; as a result, their

*O*(1) contribution is hereon ignored.

*Stationary point of type II ^{−}.* With reference to figure 2, it is clear that, depending on

**°, integrals**

*x**J*

_{1}and

*J*

_{2}may feature

*multiple*stationary points of type II

^{−}according to the second of (3.11). For this class of critical points, it is shown in electronic supplementary material, appendix A(a) that

*r*

_{1/2}and their bounds are detailed in the electronic supplementary material, appendix A(b) (see for instance, figure S15). From (3.18), it is seen that the non-uniform asymptotic expansion (3.8) breaks down as

*r*→

*r*

_{1/2}, which in physical terms corresponds to

**° straddling a**

*x**caustic region*[33]. Electronic supplementary material, appendix A(c) demonstrates that in this case the corank of

*A*

_{pq}approaches either 1 or 2, depending on

**relative to the orthonormal basis (**

*d*

*a*_{1},

*a*_{2},

**)—given by the principal directions and the outward normal to**

*n**S*

^{f}at

*ζ*^{−}

_{II}. On denoting by

*r*=

*r*

_{1/2}where (3.8) fails, the ‘minus’ counterpart of (3.17) can be shown to read

*ζ*^{−}

_{II}does not interact with its neighbours in the sense that nominally

*ζ*^{−}

_{II}, and the germane interaction must be accounted for via uniform asymptotic expansion of (3.7) that is examined next.

### (c) Uniform topological sensitivity approximation in the caustic region

To frame the above discussion in a formal setting, recall that for a given i.e. fixed obstacle shape, the *bifurcation set* [34] of the phase function
*ϕ* whose distance vanishes as (** d**,

**°)→**

*x**B*

_{ϕ}.

### Lemma 3.2

*For the problem under consideration*,
*where the loci* II^{−} *and affiliated caustic distances r*_{1/2} *are specified, respectively, in* (3.13) *and electronic supplementary material, appendix A*(*a*).

### Proof.

The claim is a direct consequence of (i) definition (3.21); (ii) the completeness of the set of stationary points given by (3.11), and (iii) the fact that the only loci in (3.13) which permit singular Hessian of the phase function are those of type II^{−}. ▪

Following [35], the interaction between stationary points should be considered as soon as their diminishing distance reaches *O*(*k*^{−1/2}) (a more precise condition will be established later). Hence when, given ** d**, the sampling point approaches the bifurcation set, i.e. straddles the caustic region, the phase function is characterized by at least two

*interacting*stationary points whose analysis warrants a uniform asymptotic treatment. In the context of (3.19), this neighbourhood of interaction, as measured along ray II

^{−}, is denoted by

#### (i) Elements of the catastrophe theory

The fundamental framework for the analysis of interacting (or coalescing) stationary points is provided by the catastrophe theory [36,37], which is rooted in the notion of structurally stable bifurcations [34]. To facilitate the discussion, assume without loss of generality that the phase function has a critical point at *η*^{1}=*η*^{2}=0 so that ∇ *ϕ*|_{0}=**0**. In this setting, the theory originates from the Morse Lemma and Splitting Lemma (e.g. [38]), which guarantee the existence of a local diffeomorphism (*η*^{1},*η*^{2})→(*ϑ*^{1},*ϑ*^{2}) in a neighbourhood of the critical point such that
*ϕ*_{°} is a constant, *ψ* is a smooth function whose value and derivatives up to order two all vanish at the origin. The basic question regarding (3.23), whose first phase representation signifies the non-degenerate case examined in §3b, deals with the order of degeneracy carried by function *ψ*. This issue is resolved via the concept of *codimension*, cod(*ϕ*)=cod(*ψ*), of the phase function that can be introduced as follows. Consider first the so-called Jacobian ideal of *ϕ*, given by Δ(*ϕ*)=*g*_{1}∂*ϕ*/∂*ϑ*_{1}+*g*_{2}∂*ϕ*/∂*ϑ*_{2} for arbitrary smooth functions *g*_{1/2}, and its formal Taylor series, *j*Δ(*ϕ*). With such definitions, the codimension of *ϕ* (assuming it is finite) can be written as
*H*_{2} is the space of all power series *j*Δ(*ϕ*) is expressible in terms of monomials, cod(*ϕ*) is simply the number of missing monomials relative to those in *H*_{2}. As examined in [34], the geometric implication of (3.24) is that a small perturbation of *ϕ* with codimension *n* can produce at most *n*+1 critical points in a neighbourhood of *η*^{1}=*η*^{2}=0.

Perhaps the most powerful result of the catastrophe theory is that of *universal unfolding*, which encapsulates feasible perturbations of *ϕ* (assuming structural stability) and provides for a uniform asymptotic treatment of diffraction catastrophes in a neighbourhood of the bifurcation set. For a phase function *c*_{m} are the control parameters that vanish on *B*_{ϕ}, and *h*_{m} form a basis for *H*_{2} modulo *ψ*, since *ϕ* and *ψ* by definition share the codimension and universal unfolding.

*Diffraction scaling.* On denoting by ** Ψ**(

*c*

_{1},…,

*c*

_{M}) the canonic Fourier integral with

*k*=1 and prototypical unfolding (3.25) of the phase function, the leading-order contribution of cognate critical point to (3.7) in the neighbourhood of

*B*

_{ϕ}when

*k*≫1 can be computed (up to an

*O*(1) multiplier) by way of

*diffraction scaling*[35] as

*k*

^{μ}

**(**

*Ψ**k*

^{σ1}

*c*

_{1},…,

*k*

^{σM}

*c*

_{M}), where

*μ*is the so-called singularity index signifying the intensity of a caustic,

*c*

_{m}are

*k*-independent, and

*σ*

_{m}>0 are the measures of fringe spacings in the control directions

*c*

_{m}(see the electronic supplementary material, appendix B for details).

#### (ii) Asymptotic order of the uniform approximation

With reference to (3.23)–(3.25), table 2 provides the complete list of elementary diffraction catastrophes with cod(*ϕ*)<4 according to Thom's classification theorem [34,35], including the respective universal unfoldings (where (*ϑ*^{1},*ϑ*^{2}) are replaced by (*s*,*t*)) and diffraction scaling parameters. Note that the diffraction catastrophes with cod(*ϕ*)>3 have not been fully analysed due to their complexity [35]. To aid the high-frequency evaluation of TS, electronic supplementary material, appendix B outlines the uniform asymptotic expansion of two-dimensional Fourier integral (3.7) for each featured type of diffraction catastrophe. The main result of this summary, listed in the last column of table 2, is the (fractional) asymptotic order of the uniform expansion when applied to the TS formula (3.3). As a point of reference, one may recall that the non-uniform approximations of type II are *O*(*k*), while those of type *I* are *O*(1).

*Global shape of a scatterer.* Assuming structural stability, the type of catastrophe affiliated with given stationary point ** d**,

**°)→**

*x**B*

_{ϕ}depends on the local behaviour of the phase function, and thus on the geometry of

*S*

^{f}, in a neighbourhood of

*second-order*properties of

*S*

^{f}(synthesized via the second fundamental form) at

*third- and higher order*surface properties of

*S*=∂

*D*[35]. The principal result of this paper in terms of theorems 4.10 and 4.12, however, applies

*regardless*of this caveat—as long as the diffraction catastrophes affiliated with

*S*do not exceed three in terms of their codimension.

### (d) Topological sensitivity approximation in the neighbourhood of *S*^{f}

To complete the analysis, consider the case *S*^{f} given by (3.5). It is apparent from figure 2 that as ** x**°→

*S*

^{f}from the outside (resp. inside) there exist at least two stationary points,

*ζ*^{+}

_{II}and

*ζ*^{−}

_{I}(resp.

*ζ*^{−}

_{II}and

*ζ*^{+}

_{I}), that merge at the normal projection of

**° onto**

*x**S*

^{f}, denoted by

*x*^{⋆}. Further when

**°∈**

*x**S*

^{f}, the phase function in (3.6) assumes locally conical shape and becomes non-differentiable at

*r*=0, i.e.

**=**

*ζ***°=**

*x*

*x*^{⋆}, which is also the point where the non-exponential factors of integrands in

*J*

_{1}and

*J*

_{2}become singular. Under such circumstances, the asymptotic approximations developed in §3b,c break down, i.e. cease to represent the contribution of stationary points located in the vicinity of

*x*^{⋆}. The purpose of this section is accordingly twofold, namely to (i) identify the length scale

*ϵ*in (3.5) which preserves the validity of previously developed approximations and (ii) expose the asymptotic contribution of

*x*^{⋆}∈

*S*

^{f}to the TS (3.3) when

#### (i) Extent of N ϵ

With reference to figure 3*a*, consider without loss of generality the situation where
*x*^{⋆}∈*S*^{f} and small |ℓ|, and let ** ζ***=

*ζ*^{−}

_{II}denote the germane stationary point of type II. Next, recall the two-term extension [33] of the non-uniform approximation (3.8) which reads

*φ*(

**), where**

*ζ**A*

_{pq}=∂

^{2}

*φ*/(∂

*η*

^{p}∂

*η*

^{q});

*g*

_{ij}=(

*i*!

*j*!)

^{−1}∂

^{i+j}

*g*(

**)/(∂**

*ζ**x*

^{i}∂

*y*

^{j})|

_{ζ=ζ*}for

*g*=

*φ*,

*f*and (

*x*,

*y*) are obtained by a local diffeomorphism from (

*η*

^{1},

*η*

^{2}) so that ∂

^{2}

*φ*/∂

*x*∂

*y*=0 at

**=**

*ζ****.**

*ζ*In the context of (3.26), the idea behind exposing the characteristic length *ϵ* in (3.5) is to find a threshold value of |ℓ| beyond which |*k*^{−1}*f*_{1}/*f*_{0}|=*o*(1). For brevity of exposition, the attention is hereon focused on applying (3.26) to the component of *J*_{1} in (3.6) with phase function ** ζ**⋅

**−**

*d**r*, noting that the analysis of the remaining integrals in (3.6) yields the same result when

**°=**

*x*

*x*^{⋆}+ℓ

**(**

*n*

*x*^{⋆}). To commence the analysis, let |ℓ|=|

**°−**

*x*

*x*^{⋆}|=

*O*(

*k*

^{−1}), and let

*x*^{⋆}be located away from ∂

*S*

^{f}so that

*** in figure 3**

*ζ**a*is contained within a ball

*O*(

*k*

^{−1}) centered at

*x*^{⋆}. In the high-frequency regime, one has

*ρ*

_{1}≥

*ρ*

_{2}≫

*k*

^{−1}, where

*ρ*

_{1}and

*ρ*

_{2}are the principal radii of curvature of

*S*

^{f}at

*x*^{⋆}. As a result,

*S*

^{f}can be locally approximated (within

*Π*

_{x⋆}, drawn at

*x*^{⋆}. As shown in §3d(ii), this treatment induces

*O*(

*k*

^{−2}) error in the integration procedure.

To aid the application of (3.26), let the normal projection of ** ζ**∈

*S*

^{f}on

*Π*

_{x⋆}be specified in terms of Cartesian coordinates (

*x*,

*y*) such that: (i)

*x*^{⋆}is identified with the origin (0,0) and (ii)

*x*is parallel to the tangential component of

**, namely**

*d*

*d*_{t}=

**+|**

*d***⋅**

*d***|**

*n***(**

*n*

*x*^{⋆}). In this setting, the phase function can be approximated locally as

*k*. On computing the projection of the stationary point

*** onto**

*ζ**Π*

_{x⋆}as (

*x**,

*y**)=(−ℓ|

*d*_{t}|/

*d*

_{n},0), where

*d*

_{n}=|

**⋅**

*d***(**

*n*

*x*^{⋆})|, the reduced phase function (3.28) can be expanded about (

*x**,

*y**) in Taylor series up to the fourth order to evaluate the necessary derivatives in (3.27). After treating in a similar way the multiplier of

*π*/

*k*. As a result, the second-order term in (3.26) can be neglected, i.e. (3.8) holds, for normal distances to

*S*

^{f}of at least

*one wavelength*. One should bear in mind that, as the shadow region is approached when

*x*^{⋆}→∂

*S*

^{f}, i.e.

*** and**

*ζ*

*x*^{⋆}exceeds

*O*(

*k*

^{−1}) (figure 3

*a*). In this case, however, the situation is mitigated by the fact that the kernels in (3.6) are all proportional to

*d*

_{n}, which makes precise knowledge of the portal distance in this border region less relevant. Accordingly, the above threshold on |ℓ| is applied uniformly ∀

*x*^{⋆}∈

*S*

^{f}by stipulating

*ϵ*=

*O*(

*k*

^{−1})≥2

*π*/

*k*in (3.5).

#### (ii) Asymptotic expansion for x ∘ ∈ N ϵ

In situations where the sampling point ** x**° straddles the ‘near-boundary’ region (3.5) with

*ϵ*=

*O*(

*k*

^{−1})≥2

*π*/

*k*, the method of stationary phase ceases to apply for critical points close to the normal projection,

*x*^{⋆}, of

**° on**

*x**S*

^{f}. Further as

**°→**

*x*

*x*^{⋆}, the normal projection itself becomes a critical point owing to the loss of differentiability of the integrands in (3.6) there. This section is devoted to computing asymptotically the contribution of

*x*^{⋆}∈

*S*

^{f}(and its neighbourhood) to T(

**°,⋅,⋅) when**

*x*It is well known that the TS can be expressed as a bilinear form entailing two forward solutions for the reference domain, namely the incident field and the so-called adjoint field (e.g. [10]). In the context of figure 1, this guarantees that T(** x**°,⋅,⋅) is in fact analytic for

*r*→0 (i.e.

**°→**

*x**S*

^{f}) are the

*artefact*of rearranging (3.4) to cater for the method of stationary phase and can be dispensed with. Focusing on the component integral

*J*

_{1}in (3.3), one finds from (2.3) and (3.4) that

*r*=0. To analyse (3.29) when

*kr*=0 and reaches maximum (absolute) value at

*kr*≃0.66 i.e.

*r*=

*O*(

*k*

^{−1}). Thus, for sufficiently high

*k*the contribution of

*x*^{⋆}to

*J*

_{1}can be evaluated by approximating

*S*

^{f}via its tangent plane (

*Π*

_{x⋆}) as shown in figure 3

*a*.

To expose the error in computing the contribution of *x*^{⋆} to *J*_{1} via tangent-plane approximation, consider a generic normal section of *S*^{f} at *x*^{⋆}, and let *b* so that
*r*=*O*(*k*^{−1}), one has *ϱ*=*O*(*k*^{−1}) and *ε*=*O*(*k*^{−2}) under the premise of locally constant radius of curvature. Accordingly, it follows from (3.29) and (3.30) that
*x*^{⋆} to *J*_{1}, where *ε*=0 in (3.30). On adopting the polar coordinate system (*ϱ*,*θ*) centred at *x*^{⋆} so that direction *d*_{t}=** d**−|

**⋅**

*d***|**

*n***corresponds to**

*n**θ*=0, one finds that

**=**

*n***(**

*n*

*x*^{⋆}),

*θ*can be computed in terms of Bessel functions of the first kind, reducing the outer integral to a pair of Hankel transforms

*x*^{⋆}to

*J*

_{1}when

*x*^{⋆}to

*J*

_{2}. On the basis of (3.3), (3.34) and (3.35), one finds

*x*^{⋆}to the TS. It is perhaps not surprising that (3.36) shares the

*common multiplier*with stationary phase approximations (3.17) and (3.19), dependent on |

**⋅**

*d***| as well as the coefficients A and B—specified by the type of (impenetrable) vanishing perturbation according to table 1.**

*n*## 4. Imaging ability of the topological sensitivity indicator function

From (2.7), it is seen that for *bi-linear* form entailing two regular wave fields in the reference domain, namely the incident wave and the fundamental solution whose source is outside *π*/*k*, i.e. half that of the illuminating wave. In this setting, the key question is that of the conditions under which the most pronounced negative values of TS are localized in a narrow region ‘about the boundary’ [11] of an obstacle.

### (a) Single plane-wave incidence

To provide an explicit platform for the analysis, the foregoing asymptotic developments (assuming the hidden anomaly to be of Dirichlet type) can be synthesized by writing
_{M}(*m*) is the characteristic function equalling 1 for *m*∈*M* and 0 otherwise; recalling figure 4, *S*^{f} given by (3.5); ** d**, the contributions of type T

^{⋆}, T

^{c}and T

^{II+}are

*unique*due, respectively, to: (i) the uniqueness of the normal projection of

*S*

^{f}, (ii) premise that the hidden obstacle is convex with smooth boundary and (iii) geometrical grounds elaborated in the electronic supplementary material, appendix C. By contrast, T(

**°) may include the contribution of**

*x**multiple*isolated stationary points of type II

^{−}, as indicated by the summation symbol before T

^{II−}. In the context of (4.1), one should also mention that for

*O*(

*k*

^{−1}) from

*x*^{⋆}—accounted for via T

^{⋆}—is implicitly excluded when computing T

^{c}and T

^{II±}.

From (4.2), it is readily seen that the near-boundary contribution is *O*(*k*) i.e. commensurate with the non-uniform approximation, yet *subpar in order* relative to the asymptotic contribution of diffraction catastrophes summarized in table 2. Accordingly the high-frequency distribution of TS is, under the premise of single plane-wave incidence, asymptotically dominated by the caustics.

### (b) Full source aperture

To expose the imaging ability of the TS indicator function, consider the full source aperture companion of (2.7), namely
** d** of incident plane wave,

*Ω*is the unit sphere, and the dependence of T on

**is implicit.**

*d*

### Proposition 4.1

*For given* *every boundary point* ** ζ**∈

*S becomes stationary point of type*II

*for some unique incident direction*

**=**

*d****(**

*d***°,**

*x***)**

*ζ**provided that*

### Proof.

From (3.11), one finds that for stationary points of type II, ** d*** must satisfy

*ζ*^{±}

_{II}∈

*S*

^{f}(

***). A contraction of (4.4) with**

*d***(**

*n*

*ζ*^{±}

_{II}) yields

**∈**

*ζ**S*is a stationary point of type II for pair (

**°,**

*x****), and [**

*d***−2**

*I***⊗**

*n***] is an (orthogonal) reflection matrix. The uniqueness of**

*n**** is then verified by contradiction noting that**

*d***(**

*n***) is single-valued. ▪**

*ζ*

### Remark 4.2

The uniqueness of ** d*** ceases at boundary points

*ζ*_{°}∈

*S*, where

*ζ*_{°}is a stationary point of type II

^{±}when

### Remark 4.3

For ** x**°∈

*D*, one has

**∈**

*ζ**S*. As a result, the stationarity type of every boundary point

**when**

*ζ***=**

*d****(**

*d***°,**

*x***) is II**

*ζ*^{−}. When

*S*of a Dirichlet obstacle can be split into subsets

*ζ*_{°}where

***⋅**

*d***(**

*n*

*ζ*_{°})=0 (figure 5

*a*). When

### Remark 4.4

For given *B*_{ϕ}(** d**,

**°) on the unit sphere spanned by**

*x***is a union of smooth curves and points where such curves join, intersect, or terminate as indicated in figure 5**

*d**b*. As

*dim*(

*Ω*)=2, the only diffraction catastrophe affiliated with the

*curves*in

*B*

_{ϕ}is of type fold (cod(

*ϕ*)=1), while the higher order catastrophes (cod(

*ϕ*)>1) appear as

*points*[35] on

*Ω*. In general,

*B*

_{ϕ}is contained within an open neighbourhood

*B*

_{ϕ}and

*level set*|

**|=0 and**

*c**neighbourhood*|

**|<**

*c**k*

^{−υ}, where

*υ*>0 is a catastrophe-specific scaling parameter to be specified later. For completeness, it is noted that

*B*

_{ϕ}is

*closed*for the assumption to the contrary would require cod(

*ϕ*)=0 [40]. In the context of (4.5) relating (for given

**°)**

*x***=**

*d**** to the stationary point(s)**

*d***∈**

*ζ**S*, the subset of

*S*corresponding to

In the light of the above remarks, one may observe that a *discrete* set of critical points contributing to T(** x**°) in the case of a single incident wave, see (4.1), transitions in the course of full-aperture illumination into a

*continuous*set

*S*of all boundary points contributing to

^{II+}with respect to

**as it contributes to (4.3). Next, recall that (4.5) provides the map relating**

*d***=**

*d**** to the solid angle of a boundary point with respect to**

*d***°, namely**

*x**one-to-one*on account of the uniqueness of

^{⋆}, see §3d. Hence

*S*

^{+}

_{II}in computing (4.6) via the concept of van der Corput neutralizers [32,41] (see also the electronic supplementary material, appendix B). This tool is implicitly used in all cases where the partitioning of a domain of integration is in order.

The same change of variable can be applied to the integral over ** d**. In this case, however, (4.5) is not one-to-one—which signifies the multiplicity of

*ζ*^{−}

_{II}and thus inherently accounts for the summation over T

^{II−}.

### Lemma 4.5

*By the way of* (4.1), (4.3) *and relationship* *the full-aperture distribution* *can be recast as*
*where* ** d***

*solves (4.5)*;

*is the ‘full-aperture’ neighbourhood of S constructed from (3.5)*; T

^{⋆}=0

*for*

**⋅**

*d***(**

*n*

*x*^{⋆})>0;

*and*

*is a ball of radius O*(

*k*

^{−1})

*centred at the normal projection*

*x*^{⋆}

*of*

**°**

*x**on S (figure 3*T

*a*). Geometrically, the respective support of^{⋆},T

^{c}

*and*T

^{±}

*in*(4.7)

*can be described as unit hemisphere, a small neighbourhood of the bifurcation set B*

_{ϕ}

*on the unit sphere, and the boundary of the scatterer excluding its subsets contributing to*T

^{⋆}

*and*T

^{c}.

#### (i) Contribution of non-degenerate stationary points

### Proposition 4.6

*The contribution of isolated stationary points to* *in* (4.7) *scales as*
*for sufficiently large k, assuming the codimension of phase singularities in the featured integral not to exceed three*.

### Proof.

By the way of (3.17) and (3.19), the left-hand side of (4.9) can be rewritten as a Fourier integral
*S*^{±} is such that *r*>2*πk*^{−1} thanks to (4.8), and
*r*|** d***⋅

**|**

*n*^{2}which satisfy

*ρ*

_{1/2}are the principal radii of curvature of

*S*

^{±}at

**. From the definition of**

*ζ**S*

^{±}in (4.8), the relevant roots of (4.11) represent the normal projection of

**° on**

*x**S*, i.e.

*S*

^{+}, the solution of (4.12) is unique due to the convexity and smoothness of

*S*, whereby the phase function in this case possesses a single isolated stationary point. From (3.8), one accordingly finds that the integral over

*S*

^{+}in (4.9) scales as

*O*(1).

On the other hand, the normal projection of ** x**° on

*S*

^{−}is generally not unique. In this case, it can be shown via (4.8) and (4.11) that the Hessian of

*r*|

**⋅**

*d***|**

*n*^{2}becomes singular at a critical point

**∈**

*ζ**S*

^{−}solving (4.12) only if

*S*

^{−}in (4.9), on accounting for catastrophes where both (4.12) and (4.13) hold, may include contributions of orders shown in table 3. ▪

#### (ii) Contribution of diffraction catastrophes

### Proposition 4.7

*For sufficiently large k, the contribution of caustics to* *in* (4.7) *behaves as*

### Proof.

The idea behind establishing (4.14) is to expose the measure of the vanishing support, *b*. In particular as ** d** (for given

**°) leaves the neighbourhood**

*x**uniform*approximation to either its non-uniform counterpart, or zero—on the dark side of some caustics (e.g. fold) due to the absence of real stationary points [33]. On denoting

*b*

_{m}=

*k*

^{σm}

*c*

_{m}, such transition in (4.15) occurs when

*sufficient condition*for estimating the extent of

The next step in the analysis is to establish (for given ** x**°) a

*linearized*relationship between |

**| and dist(**

*c***,**

*d**B*

_{ϕ}) on the unit sphere, in a small neighbourhood of the bifurcation set. In the context of figure 5

*b*, it is recalled that the fold caustics (cod(

*ϕ*)=1) translate into smooth non-intersecting curves in

*B*

_{ϕ}⊂

*Ω*, while the catastrophes of higher codimension are projected as points in

*B*

_{ϕ}. In this setting, dist(

**,**

*d**B*

_{ϕ}) is identified as the normal spherical distance to a curve (resp. spherical distance to a point) when cod(

*ϕ*)=1 (resp. cod(

*ϕ*)>1). On writing the sought relationship as |

**|=**

*c**V*⋅dist(

**,**

*d**B*

_{ϕ}), one finds from (4.16) that

*k*, noting that

*V*>0 since the bifurcation set

*B*

_{ϕ}⊂

*Ω*is

*closed*(see remark 4.4).

To estimate *V* , consider first a *fold* bifurcation point ** d**°∈

*B*

_{ϕ}for given

**°, and let**

*x**** ∈**

*ζ**S*denote the affiliated critical point on the boundary of the scatterer. In this case

**=**

*c**c*, and the Hessian of

*ϕ*is of corank one. On account of the Splitting Lemma (3.23) and Taylor expansion of

*ϕ*(

**) about**

*ζ**** , there exist local surface coordinates (**

*ζ**σ*,

*τ*) such that

*ϕ*

_{°}=

*ϕ*(

***);**

*ζ**ϕ*

_{°}′′′ are

*O*(1), and

*σ*,

*τ*) at

***, used to describe**

*ζ***to the leading order. The objective is to find the variation in**

*ζ**c*due to infinitesimal perturbation d

**⊥**

*d***°. Using (4.18) and definition**

*d**ϕ*=

**⋅**

*ζ***−**

*d**r*where

*r*=|

**°−**

*x***|, one finds**

*ζ**fold*universal unfolding as in table 2, one finds from (4.19) via mapping

*t*=(|

*ϕ*

_{°}′′′|/2)

^{1/3}

*τ*that

**is parallel to**

*d**c*remains zero to the leading order. This shows that

*tangent*to the fold curve at

**°∈**

*d**B*

_{ϕ}. Subsequently, the width of a stripe-like region

*B*

_{ϕ}is exposed by considering d

**in the plane containing**

*d***° and**

*d***|=dist(**

*d***,**

*d**B*

_{ϕ}),

**=**

*n***(**

*n****), and**

*ζ**ϑ*is the angle between

**° and**

*d***. From (4.17), (4.20) and table 2, one finds an upper-bound estimate**

*n*^{c}scale with |

**°⋅**

*d***|, so that the situations of widening**

*n***°⋅**

*d***|→0 pose no problem in terms of the contribution of the fold catastrophes to (4.14).**

*n*From (4.20), it is seen that the sole situation precluding *V* =*O*(1) is |** d**°⋅

**(**

*n****)|≪1. As shown in the electronic supplementary material, appendix A(b), this requires that the distance between**

*ζ***° and the critical point, |**

*x***°−**

*x****|, behaves as**

*ζ**O*(|

**°⋅**

*d***|**

*n*^{±1}), see also the electronic supplementary material, figure S15. Owing to the regularity of

*S*, however, catastrophes with cod(

*ϕ*)>1 cannot occur arbitrarily close to

*S*, while the sampling points where |

**°−**

*x****|≫1 are outside of**

*ζ**V*=

*O*(1) for higher codimension catastrophes and consequently

*isolated points*(figure 5

*b*). The claim (4.14) then follows from the scaling of T

^{c}in table 2, (3.6), (4.21) and (4.22). For completeness, the effect on (4.14) due to error of the Kirchhoff approximation (3.1) is examined in the electronic supplementary material, appendix D. ▪

#### (iii) Contribution of nearby critical points for x ∘ ∈ N ˘ ϵ

### Proposition 4.8

*For* *the contribution of nearby critical points to* *in* (4.7) *is given by*
*for sufficiently large k, where* ℓ=(** x**°−

*x*^{⋆})⋅

**(**

*n*

*x*^{⋆})

*is the signed normal distance between*

**°**

*x**and the boundary of the scatterer*.

### Proof.

The single-incident-wave expression for T^{⋆}(** x**°,⋅,⋅), given by (3.36), is explicit and permits direct integration with respect to

**. Recall that T**

*d*^{⋆}=0 for

**⋅**

*d***(**

*n*

*x*^{⋆})≥0 thanks to (3.1), whereby the effective integration support in (4.23) is a hemisphere. By taking −

**(**

*n*

*x*^{⋆}) as the zenith direction of the spherical coordinate system describing

*Ω*, it is evident from (3.36) that T

^{⋆}is exclusively a function of the zenith angle

*k*ℓ so that

### Remark 4.9

As discussed in remark 3.1, proposition 4.8 is established under the premise that a hidden obstacle *D* is sound-soft. Accordingly, it is of interest to examine the near-boundary variation (4.23) when (A,B)=(0,1), i.e. when the vanishing perturbation in (2.6) is likewise of Dirichlet type (table 1). This behaviour is shown in figure 6*a*, which plots (4.23) with (A,B)=(0,1) versus *k*ℓ in a neighbourhood of ∂*D*. As can be seen from the graph, the leading contribution of T^{⋆} to *D* and (ii) attains extreme negative (resp. positive) value at its first peak inside (resp. outside) the obstacle, at a normal distance of |*k*ℓ|<*π*/2 from the boundary. For completeness, figure 6*b* plots the corresponding distribution of (4.23) assuming *a*. This observation will motivate the proposed TS algorithm for high-frequency obstacle reconstruction.

The foregoing developments are now concluded with the main result of this work.

### Theorem 4.10

*Consider the inverse scattering problem for a convex Dirichlet obstacle D as in figure 1 with far-field sensory data. For sufficiently large k, the full-source-aperture distribution of the scaled TS indicator* (2.6) *behaves as*
*under the premise of diffraction catastrophes with codimension less than four, where* *is a* 2*ϵ-thick shell* (*for some ϵ*=*O*(*k*^{−1})≥2π/*k*) *with mid-plane ∂D*; ℓ=(** x**°

*−*

*x*^{⋆})⋅

**(**

*n*

*x*^{⋆})

*is the signed normal distance between*

**°**

*x**and ∂D, and the coefficient pair*(A,B)

*takes values as in table 1 depending of the type of*(

*impenetrable*)

*vanishing perturbation used to probe the domain.*

### Proof.

The claim is a direct consequence of lemma 4.5 and propositions 4.6–4.8. As shown in the electronic supplementary material, appendix D, the result is resilient to the error due to Kirchhoff approximation (3.1). ▪

### Remark 4.11

Note that (4.25) ensures the high-frequency reconstruction of a Dirichlet obstacle assuming full source aperture (** d**∈

*Ω*). At a glance, such requirement may appear excessive in the light of the uniqueness-of-reconstruction result for sound-soft obstacles (corollary 5.3 in [43]) with only a single incident plane wave. As examined in [44], however, the latter claim holds only for scatterers contained within a ball of radius

*R*

_{°}≃4.49/

*k*—which inherently precludes the high-wavenumber case considered herein. Indeed, the numerical results in §5 (see for instance, figure 9) demonstrate that, at wavelengths that are small relative to the obstacle size, the TS reconstruction with a single incident plane wave provides no information about the ‘dark side’ of the obstacle. On the other hand, from theorem 5.1 in [43] it follows that the above uniqueness result does hold for any wavenumber

*k*provided that the obstacle is illuminated by an infinite number of incident plane waves—which is consistent with the claim of theorem 4.10.

### (c) Neumann obstacle

For a sound-hard obstacle, the physical optics approximation [32] reads
*u*^{i}=e^{−ikx⋅d} to (2.8), followed by the use of (2.13) and (2.14) to address the component integrals over *Γ*^{obs}, results in a TS formula for Neumann obstacle that is structurally similar to (3.3). In particular, the kernel in the ‘sound-hard’ counterpart of (3.3) can be shown to (i) feature the identical phase function ** ζ**⋅

**±**

*d**r*, (ii) remain regular as

**°→**

*x**S*

^{f}and (iii) vanish on ∂

*S*

^{f}. The end result of the analysis is given by the following statement.

### Theorem 4.12

*Consider the inverse scattering problem for a convex Neumann obstacle D as in figure 1 with far-field sensory data. For sufficiently large k, the full-source-aperture distribution of the scaled TS indicator* (2.6) *behaves as**under the premise of diffraction catastrophes with codimension less than four, where* *is a* 2*ϵ-thick shell* (*for some ϵ*=*O*(*k*^{−1})≥2π/k) *with mid-plane ∂D*; ℓ=(** x**°

*−*

*x*^{⋆})⋅

**(**

*n*

*x*^{⋆})

*is the signed normal distance between*

**°**

*x**and ∂D, and the coefficient pair*(

*A*

*,*

*B*)

*takes values as in table*1

*depending of the type of (impenetrable) vanishing perturbation used to probe the domain.*

### (d) Unscaled topological sensitivity distribution

To establish a direct link with the results of earlier TS studies on impenetrable obstacles (e.g. [27]), it is of interest to rewrite (4.25) and (4.27) in terms of the original (i.e. unscaled) TS formula (2.5), and to further assume *prior knowledge* of obstacle type by letting (A,B)=(0,1) (resp. *T*(** x**°)=

*k*

^{−γ}T(

**°) according to (2.6) and table 1, one consequently finds from theorems 4.10 and 4.12 that**

*x**Dirichlet*obstacles probed by sound-soft perturbations, and

*Neumann*anomalies probed by sound-hard perturbations. As can be seen from (4.28) and (4.29), the two (unscaled) asymptotic behaviours are quite distinct, yet both localized near ∂

*D*. This can be verified by noting that the Dirichlet variation (4.28) is, up to the factor

*k*

^{−2}, given by figure 6

*a*, while the Neumann variation (4.29) is—up to the sign factor—shown in figure 6

*b*.

### (e) Reconstruction scheme

A comparison between (4.25) and (4.27) reveals that for (A,B) *fixed*, the leading-order behaviour of *changes sign* when (2.6) is applied towards the reconstruction of a *hidden Dirichlet* versus *hidden Neumann* obstacle. Accordingly, the counterpart of figure 6 for a hidden Neumann anomaly is obtained via reflection of the featured diagrams about the *k*ℓ-axis. This opens two distinct avenues towards the high-frequency TS reconstruction of impenetrable obstacles:

### Algorithm 4.13

*When the nature of D is known beforehand*, (i) *compute* *with commensurate trial parameters* ((A,B)=(0,1) *for Dirichlet anomaly*, *for Neumann anomaly*) *and* (ii) *reconstruct* ∂*D as the zero level set of* *separating its extreme negative and extreme positive values. The extreme* *values inside the reconstruction are always negative*.

### Algorithm 4.14

*When the nature of D is unknown, compute* *with* (A,B)=(0,1) *and reconstruct ∂D as the zero level set of* *separating its extreme negative and extreme positive values. When the extreme* *values to the inside of the reconstruction are negative* (*resp. positive*), *the impenetrable obstacle is of Dirichlet* (*resp. Neumann*) *type*.

In this setting, algorithm 4.14 is generally preferred due the to the facts that: (i) the obstacle type is *revealed* rather than required as prior information, and (ii) setting (A,B)=(0,1) (as opposed to *stronger localization* of the extreme *D* (figure 6).

## 5. Numerical results

A computational experiment is devised to illustrate the performance of the TS as an imaging tool in the high-frequency regime. The focus is on elucidating how does (3.3) relate to the boundary of a convex Dirichlet obstacle, and how is the numerical distribution thereof approximated by the closed-form expression (4.25) in the case of the full source aperture. The sensing arrangement is shown in figure 7*a*, where *D* is an ellipsoidal anomaly with semi-axes (0.2,0.08,0.8). In what follows, the TS distribution is computed in the obstacle's mid-section *Π* perpendicular to its major axis, assuming incident plane waves with *k*=300 (wavelength λ=0.021) propagating in direction ** d**∥

*Π*. In this case, the computation of TS is facilitated by two critical observations: (i) for

**° within the square ‘patch’ shown in figure 7**

*x**a*, the critical points on

*S*

^{f}are confined to

*S*

^{f}∩

*Π*and (ii) the germane catastrophes are of either type fold or cusp, i.e. cod(

*ϕ*)≤2. Note, however, that the caustics of higher codimension may occur in out-of-plane situations—when either

**°∉**

*x**Π*—see figure 7

*b*for an example projection of the bifurcation set

*B*

_{ϕ}(

**,**

*d***°) on**

*x**Ω*.

### (a) Applicability of Kirchhoff approximation

One may recall that the implicit assumptions behind (3.3)—used hereon as the basis for numerical evaluation of TS—are that: (i) the sensory data are of far-field type (see remark 2.1) and (ii) the Kirchhoff approximation (3.1) applies. In this setting, it is of interest to assess the accuracy of (3.1) for the testing configuration described above. Owing to lack of suitable three-dimensional solutions, the validation is performed in a two-dimensional setting, assuming scattering by a sound-soft *cylinder* whose cross-section equals *S*∩*Π* in figure 7. With such premise, figure 8 compares the analytical solution [45] for the far-field pattern,

### (b) Single plane-wave incidence

With the aid of the high-frequency approximations described in §3 and the electronic supplementary material, appendix B, the TS field (4.1) is computed via the following steps: (i) the 0.5×0.5 square computational domain within *Π* (figure 7*a*) is discretized by 10^{6} pixels, nearly 42 per wavelength; (ii) the boundary curve *S*^{f}(** d**)∩

*Π*is split into 10

^{4}segments centred at

*ζ*^{1}, the contribution of

*ζ*^{n}to (4.1) is computed (via either near-boundary, uniform or non-uniform approximation) along rays II

^{±}∈

*Π*, and accordingly used to ‘paint’ the pixels. In doing so, the use is made of the van der Corput neutralizers [32] to prevent double-counting of individual contributions. Assuming the ellipsoidal anomaly to be of Dirichlet type, the resulting TS map is shown in figure 9

*a*, which clearly reflects the presence of fold- and cusp-type caustics. For completeness, figure 9

*b*plots the corresponding diagram obtained via numerical integration of (3.3), while figure 9

*c*compares the two estimates along an example ray II

^{−}. As can be seen from the panel, the near-boundary, uniform and non-uniform approximations smoothly transition into one another and overlap with the numerical solution.

In the context of remark 4.11, one may observe from figure 9 that obstacle illumination by a single incident wave yields no information about the ‘dark side’ of the anomaly thanks to the fact that the scattered field vanishes (to the leading order) there, see (3.1). Recalling Corollary 5.3 in [43], on the other hand, it is further noted (using the length scale in figure 9) that *D* would have to be contained within a ball of radius *R*_{°}≃4.49/*k*≃0.015 [44] in order to guarantee the uniqueness of obstacle reconstruction with only a single incident plane wave.

### (c) Partial and full source aperture

In what follows, the partial- and full-source-aperture simulations of TS are effected by first (i) integrating (3.4) numerically to obtain (3.3)—as a function of ** x**°—for given

**, and then (ii) integrating the latter result (also numerically) with respect to**

*d***over a prescribed subset of**

*d**Ω*. In the context of the single-incident-wave example in figure 9, it is first of interest to integrate T(

**°) with respect to**

*x***∥**

*d**Π*, i.e. with respect to the

*in-plane*angle of incidence

*θ*shown in figure 9

*a*. On denoting for brevity

^{⋆}, T

^{c}and T

^{II±}to

*O*(

*k*),

*O*(

*k*

^{α}) and

*O*(

*k*

^{μ}), where

*in-plane*aperture. Note that (i) the bright sector of the unit circle in each panel depicts the source aperture, (ii) the TS distributions are thresholded at 45%, (iii) the bottom right panel plots

**spanning the unit circle.**

*d*In practice, one may expect to obtain a satisfactory TS reconstruction—at least at lower frequencies—with only a limited number of incident plane waves (e.g. [12,28]). To examine this possibility, figure 11 plots four ‘coarse’ approximations of *N*≤32 incident wave directions (indicated on the unit circle) spanning [0,2*π*]. As can be seen from the display, the high-frequency TS reconstruction is notably sensitive to the density of plane-wave illumination, primarily due to the apparent effect of *caustics*. This finding appears to be consistent with the high-frequency numerical study in [8] on the TS reconstruction of Neumann obstacles in

For completeness, the reconstruction of a Dirichlet obstacle by *a*,*d*), (*b*,e) and (*d*,*f*) plot, respectively, ^{⋆} only. Note that the featured images are obtained by adopting *Algorithm 1*, which samples each anomaly with physically compatible vanishing obstacle. As can be seen from figure 12e, this leads to an apparent ‘smearing’ in the case of a Neumann obstacle. By contrast, its image obtained via *Algorithm 2*, i.e. using *b*—and thus better localized.

To provide the *full-source-aperture* counterpart of the result in figure 12*c*—computed at boundary point *x*^{⋆}=(0.178,0.036,0), figure 13 compares the analytical expression (4.25) with a numerical estimate of *πN*^{−1}×(3.3) for *N*=512 incident plane-wave directions, uniformly distributed over *Ω*. For generality, the comparison is made at both in-plane boundary point *x*^{⋆}=(0.178,0.036,0) (*a*), and its out-of-plane companion *x*^{⋆}=(0.114,0.046,0.470) (*b*). Irrespective of the boundary point, the numerical result closely follows (4.25) at both *k*=300 and *k*=600, showing visibly better agreement in the latter case.

To conclude the study, algorithm 4.13 is applied to identify the boundary of a circular hole (Neumann obstacle) in an aluminium plate from the recent set of *elastodynamic* experiments [21]. In this case, elastic waves are propagated in a *bounded* domain shown in figure 14*a* and monitored along its top and side edges. The incident waves are generated by a piezoelectric transducer, placed sequentially at five locations indicated in the panel, such that the ratio between the wavelength and the obstacle size is 0.85. Thus, the testing configuration is incompatible with this analysis in several aspects, including (i) dimensionality of the problem, (ii) type of the governing equation, (iii) geometry of the anomaly-free domain, (iv) probing wavelength, and (v) source aperture. Nonetheless, the reconstruction of a circular hole in panel (*c*), obtained by applying algorithm 4.13 to the TS distribution [21] shown in (*b*), is rather satisfactory.

## 6. Conclusion

In this work, it is shown why the TS may work as a non-iterative tool for the waveform tomography of finite-sized anomalies in the short wavelength regime. The analysis confirms previous numerical and experimental findings to this effect, which have so far eluded rigorous justification. To establish the claim, it is assumed that anomaly is convex and impenetrable, and that the sensory data are of the far-field type. Making use of the multipole expansion and Kirchhoff approximation, the TS indicator function is first expressed as a surface Fourier integral over the illuminated part of obstacle's boundary. Under the high-wavenumber hypothesis, the latter is pruned to three asymptotic essentials, namely (i) the near-boundary approximation for sampling points within few wavelengths from the illuminated surface of an anomaly, (ii) diffraction catastrophes (of codimension <4) for sampling points near caustic surfaces, lines and points and (iii) stationary phase approximation in the remainder of the sampled region. In the case of the full source aperture, it is shown via catastrophe theory that the TS is asymptotically dominated by the explicit near-boundary term. This unveils the new reconstruction logic at short wavelengths, where the anomaly's boundary is obtained as a zero level set of the TS field separating its extreme negative and extreme positive values while its character—if unknown beforehand—is exposed from the sign of the near-boundary variation. The analysis inherently lends itself to the treatment of diffraction catastrophes with higher codimension (≥4) for which uniform approximation may become available. However, extensions of the study to penetrable, non-convex or multiple scatterers remain open question due to lack of explicit (Kirchhoff-type) approximations for the scattered field on the boundary of such obstacles.

## Competing interests

We declare we have no competing interests.

## Funding

The support provided by the US Department of Energy via NEUP grant no. 10-862 and the University of Minnesota Supercomputing institute is kindly acknowledged.

## Acknowledgements

Special thanks are due to Sir Michael Berry and David Farrelly for their input during the course of this investigation.

- Received March 16, 2015.
- Accepted May 15, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.