Characterisation of Multiple Conducting Permeable Objects in Metal Detection by Polarizability Tensors

Realistic applications in metal detection involve multiple inhomogeneous conducting permeable objects and the aim of this paper is to characterise such objects by polarizability tensors. We show that, for the eddy current model, the leading order terms for the perturbation in the magnetic field, due to the presence of $N$ small conducting permeable homogeneous inclusions, comprises of a sum of $N$ terms with each containing a complex symmetric rank 2 polarizability tensor. Each tensor contains information about the shape and material properties of one of the objects and is independent of its position. The asymptotic expansion we obtain extends a previously known result for a single isolated object and applies in situations where the object sizes are small and the objects are sufficiently well separated. We also obtain a second expansion that describes the perturbed magnetic field for inhomogeneous and closely spaced objects, which again characterises the objects by a complex symmetric rank 2 tensor. The tensor's coefficients can be computed by solving a vector valued transmission problem and we include numerical examples to illustrate the agreement between the asymptotic formula describing the perturbed fields and the numerical prediction. We also include algorithms for the localisation and identification of multiple inhomogeneous objects.


Introduction
There is considerable interest in being able to locate and characterise multiple conducting permeable objects from measurements of mutual inductance between a transmitting and a receiving coil, where the coupling is inductive. An obvious example is in metal detection where the goal is to identify and locate the multiple objects present in a low conducting background. Applications include security screening, archaeological digs, ensuring food safety as well as the search for landmines and unexploded ordnance and landmines. Other applications include magnetic induction tomography for medical imaging applications and monitoring of corrosion of steel reinforcement in concrete structures. In all these practical applications, the need to locate and distinguish between multiple conducting permeable inclusions is common. This includes benign situations, such as coins and keys accidentally left in a pocket during a security search or a treasure hunter becoming lucky and discovering a hoard of Roman coins, as well as threat situations, where the risks need to be clearly identified from the background clutter. For example, in the case of searching for unexploded landmines, the ground can be contaminated by ring-pulls, coins and other metallic shrapnel, which makes the process of clearing them very slow as each metallic object needs to be dug up with care. Furthermore, conducting objects are also often inhomogeneous and made up of several different metals. For instance, the barrel of a gun is invariably steel while the frame could be a lighter alloy, jacketed bullets have a lead shot and a brass jacket and modern coins often consist of a cheaper metal encased in nickel or brass alloy. Thus, in practical metal detection applications, it is important to be able to describe both multiple objects and inhomogeneous objects. Magnetic polarizability tensors (MPTs) hold considerable promise for the low-cost characterisation in metal detection. An asymptotic expansion describing the perturbed magnetic field due to the presence of a small conducting permeable object has been obtained by Ammari, Chen, Chen, Garnier and Volkov [2], which characterises the object in terms of a rank 4 tensor. Ledger and Lionheart have shown that this asymptotic expansion simplifies for orthonormal coordinates and allows a conducting permeable object to be characterised by a complex symmetric rank 2 MPT with an explicit expression for its 6 coefficients [10]. Ledger and Lionheart have also investigated the properties of this tensor [11] and they have written the article [12] to explain these developments to the electrical engineering community as well as to show how it applies in several realistic situations. In [13] they have obtained a complete asymptotic expansion of the magnetic field, which characterises the object in terms of a new class of generalised magnetic polarizability tensors (GMPTs), the rank 2 MPT being the simplest case. The availability of an explicit formula for the MPT's coefficients, and its improved understanding, allows new algorithms for object location and identification to be designed e.g. [3]. Electrical engineers have applied MPTs to a range of practical metal detection applications, including walk through metal detectors, in line scanners and demining e.g. [8,18,17,4,16,7,20], see also our article [12] for a recent review, but without knowledge of the explicit formula described above. Engineers have made a prediction of the form of the response for multiple objects e.g. [5], but without an explicit criteria on the size or the distance between the objects in order for the approximation to hold. Grzegorczyk, Barrowes, Shubitidze, Fernández and O'Neill have applied a time domain approach to classify multiple unexploded ordinance using descriptions related to MPTs [9]. Davidson, Abel-Rehim, Hu, Marsh, O'Toole and Peyton have made measurements of MPTs for inhomogeneous US coins [6] and Yin, Li, Withers and Peyton have also made measurements to characterise inhomogeneous aluminium/carbon-fibre reinforced plastic sheets [19]. But, in all cases, without an explicit formula. Our work has the following novelties: Firstly, we characterise rigidly joined collections of different metals (i.e. metals touching or held in that configuration by a non-conducting material) by MPTs overcoming a deficiency of our previous work. Secondly, we find that the frequency spectra of the eigenvalues of the real and imaginary parts of the MPT for an inhomogeneous object exhibit multiple non-stationary inflection points and maxima, respectively, and the number of these gives an upper bound on the number of materials making up the object. To achieve this, we revisit the asymptotic formula of Ammari et al. [2] and our previous work [10] and extend it to treat multiple objects by describing the perturbed magnetic field as a sum of terms involving the MPTs associated with each of the inclusions. We also provide a criteria based on the distance between the objects, which determines the situations in which the expression will hold. We derive a second asymptotic expansion that describes the perturbed magnetic field in the case of inhomogeneous objects and, as a corollary, this also describes the magnetic field perturbation in the case of closely spaced objects. In each case, we provide new explicit formulae for the MPTs. We also present algorithms for the localisation and characterisation of objects, which extends those for the isolated object case [2]. The paper is organised as follows: In Section 2, the characterisation of a single homogeneous object is briefly reviewed. Section 3 presents our main results for characterising multiple and inhomogeneous objects by MPTs. Sections 4 and 5 contain the details of the proofs for our main results. In Section 6, we present numerical results to demonstrate the accuracy of the asymptotic formulae and presents results of algorithms for the localisation and identification of multiple (inhomogeneous) objects.

Characterisation of a Single Conducting Permeable Object
We begin by recalling known results for the characterisation of a single homogenous conducting permeable object. Following [2,10] we describe a single inclusion by B α :" αB`z, which means that it can be thought of a unit-sized object B located at the origin, scaled by α and translated by z. We assume the background is non-conducting and non-permeable and introduce the position dependent conductivity and permeability as where µ 0 is the permeability of free space and B c α :" R 3 zB α . For metal detection, the relevant mathematical model is the eddy current approximation of Maxwell's equations since σ˚is large and the angular frequency ω " 2πf is small (see Ammari, Buffa and Nédélec [1] for a detailed justification). Here the electric and magnetic interaction fields, E α and H α , respectively, satisfy the curl equations in R 3 and decay as Op|x|´1q for |x| Ñ 8. In the above J 0 is an external current source with support in B c α . In absence of an object, the background fields E 0 and H 0 satisfy (1) with α " 0. The task is to find an economical description for the perturbed magnetic field pH α´H 0 qpxq due to the presence of B α , which characterises the object's shape and material parameters by a small number of parameters separately to its location z. For x away from B α , the leading order term in an asymptotic expansion for pH α´H 0 qpxq as α Ñ 0 has been derived by Ammari et al. [2]. We have shown that this reduces to the simpler form [10, 12] 1 pH α´H 0 qpxq i "pD 2 x Gpx, zqq ij pMrαBsq jk pH 0 pzqq k`p Rpxqq i " 1 4πr 3 p3r br´Iq ij pMrαBsq jk pH 0 pzqq k`p Rpxqq i .
In the above, Gpx, zq :" 1{p4π|x´z|q is the free space Laplace Green's function, r :" x´z, r " |r| andr " r{r and I is the rank 2 identity tensor. The term Rpxq quantifies the remainder and it is known that |R| ď Cα 4 }H 0 } W 2,8 pBαq . The result holds when ν :" σ˚µ 0 ωα 2 " Op1q (this includes the case of fixed σ˚, µ˚, ω as α Ñ 0) and requires that the background field be analytic in the object. Note that (2) involves the evaluation of the background field within the object, (usually at it's centre) i.e. H 0 pzq. The complex symmetric rank 2 tensor MrαBs :"´CrαBs`N rαBs in this asymptotic expansion, which depends on ω, σ˚, µ˚{µ 0 , α and the shape of B, but is independent of z, is the MPT, and its coefficients can be computed from pN rαBsq jk :"α 3ˆ1´µ 0 µ˚˙ż Bˆe These, in turn, rely on the vectoral solutions θ k , k " 1, 2, 3, to the transmission problem where r¨s Γ denotes the jump of the function over Γ and ξ is measured from an origin chosen to be in B. In [11] we have presented numerical results for the frequency behaviour of the coefficients of M for a variety of simply and multiply connected objects. These have been obtained by applying a hp-finite element method to solve (4) for θ k , k " 1, 2, 3, and then to compute M using (3). Our previously presented results have exhibited excellent agreement with for MPTs previously presented in the electrical engineering literature. Pratical applications of the asymptotic formula have been discussed in [12].

Main Results
The asymptotic formula given in (2) is for a single isolated homogenous object. But, as described in the introduction, for realistic metal detection scenarios, measurements of the perturbed magnetic field often relate to field changes caused by the presence of multiple objects or inhomogeneous objects. The purpose of this work is to extend the description to the cases of well separated multiple homogeneous objects and inhomogeneous objects. As a result of corollary, our second main result also describes the case of objects that of objects that are closely spaced. Below we summarise the main results of our paper.

Multiple Homogeneous Objects that are Sufficiently Well Spaced
We consider N homogenous conducting permeable objects of the form pB α q pnq " α pnq B pnqz pnq 2 with Lipschitz boundaries where, for the nth object, B pnq denotes a corresponding unit sized object located at the origin, α pnq denotes the object's size and z pnq the object's translation from the origin. The union of all objects is B α :" Ť N n"1 pB α q pnq where we use a bold subscript α to denote the possibility that each object in the collection can have a different size. We also employ the same notation for the fields E α and H α , which satisfy (1). An illustration of a typical configuration is shown in Figure 1. In this figure, there are N " 3 objects, which are the spheres pB α q pnq " α pnq B pnq`zpnq , n " 1, 2, 3, where, for the nth object, α pnq is its size (here its radius) and z pnq is its translation from the origin. In the presented case B " B p1q " B p2q " B p3q is a unit sphere located at the origin although, in practice, the objects do not need to have the same shape.
We introduce ν min ď ν pnq :" ωµ 0 σ pnq pα pnq q 2 ď ν max and set α min " min n"1,...,N α pnq , α max " max n"1,...,N α pnq and require that the parameters of the inclusions be such that ν max " Op1q, 2 Note no summation over n is implied. Figure 1: Illustration of a typical situation of N " 3 objects with B α " Ť N n"1 pB α q pnq " α pnq B pnq`zpnq such that they are not closely spaced where each object pB α q pnq is a sphere, α pnq is the radius of the nth sphere, z pnq describes the translation of the nth sphere from the origin and B " B p1q " B p2q " B p3q is a unit sphere positioned at the origin. which implies that ν pnq " Op1q. The task is then to provide a low-cost description of pH α´H 0 qpxq for x away from B α . This is accomplished through the following result.
Proof. The result follows from by using a tensor representation of the asymptotic formula in Theorem 4.8, which is an extension of Theorem 3.2 obtained in [2] for N sufficiently well spaced objects. A tensor representation of this result leads to each of the N objects being characterised by a rank 4 tensor. Then, by considering each object in turn and repeating the same arguments as in Theorem 3.1 in [10], which exploits the skew symmetries of the tensor coefficients, the result stated in (5)

Single Inhomogeneous Object
In this case, B α :" Ť N n"1 B pnq α " α Ť N n"1 B pnq`z " αB`z describes a single object comprised of N constitute parts, B pnq α , such that there is a single common size parameter α, the configuration B contains the origin, and z is a single translation, as illustrated in Figure 2. Notice that for the inhomogeneous case we use B pnq α rather than pB α q pnq as α is the same for all n and we revert to the use of non-bold α subscripts for the fields E α and H α , which satisfy (1). The material parameters of the constitute parts of the object B α are where B c α :" R 3 zB α and we drop the subscript α on µ and σ when considering the object B. We redefine ν min ď ν pnq :" ωµσ pnq α 2 ď ν max with the same requirements on ν max as before. The task is then to provide a low-cost description of pH α´H 0 qpxq for x away from B α . This is accomplished through the following result. Theorem 3.3. For an inhomogeneous object B α " αB`z made up of N constitute parts with parameters such that ν min ď ν pnq ď ν max the perturbed magnetic field at positions x away from B α satisfies pH α´H 0 qpxq i "pD 2 x Gpx, zqq ij pM rαBsq jk pH 0 pzqq k`p Rpxqq i , where  which, in turn, rely on the vectoral solutions θ k , k " 1, 2, 3, to the transmission problem where ξ is measured from the centre of B and, in this case, Γ :" BB Y tBB pnq X BB pnq , n, m " 1, . . . , N, n ‰ mu.
Proof. The result follows from by using a tensor representation of the asymptotic formula in Theorem 5.3, which is an extension of Theorem 3.2 obtained in [2] for an homogeneous object to the inhomogeneous case. Using a tensor representation of this result leads to the object being characterised in terms of a rank 4 tensor. Then, by repeating the same arguments as in Theorem 3.1 in [10], which exploits the skew symmetries of the tensor coefficients, the result stated in (8) is obtained. The symmetry of MrαBs follows from repeating the arguments in Lemma 4.4 in [10].
(1) Figure 3: Illustration of a typical situation of N " 3 closely spaced objects of the form B α " where each object is a sphere, α is a single scaling parameter, z describes their translation of the configuration from the origin. Corollary 3.5. Theorem 3.3 also immediately applies to objects that are closely spaced and, in this case, B α " αB`z implies a single size parameter α and a single translation z for the configuration B. An illustration of a typical configuration is shown in Figure 3. In this figure, there are N " 3 objects consisting of three spheres configured such that they scale and translate together according to α and z, respectively, and, in this case, B is the combined configuration of three (larger) spheres with different radii and with centres located away from the origin.
Remark 3.6. The applicability of Theorem 3.3 to closely spaced objects is expected to be limited since, in order to compute the characterisation, prior knowledge of the multiple object configuration (ie location and orientation with respect to each other) is required, which, in practice, will not be the case. The formula also requires that the objects be closely spaced as there is a single scaling parameter and single translation that describes the configuration, but prior knowledge of the location of the configuration is not required. Instead, this result is expected to be of more practical value in the description of inhomogeneous objects where the configuration of the different regions of an object will be known in advance.
Remark 3.7. The translation invariance of the MPT coefficients described by Proposition 5.2 in [3] and the tensor transformation rules described in the proof of Theorem 3.1 in [10] carry over immediately to the rank 2 MPTs defined in (6) and (9). 4 Results for the Proof of Theorem 3.1

Elimination of the Current Source
Recall from [2] that and the weak solution for the interaction field is : where p¨,¨q Ω denotes the standard L 2 inner product over Ω. In a departure from [2], we have, for multiple objects, that Noting that the weak solution for the background field is: we can write: Find E α PX α such that which eliminates the current source. We also obtain that

Energy Estimates
In [2] a vector field F pxq was introduced such that its curl is equal to the first two terms of a Taylor's series expansion of ∇ˆE 0 about z for |x´z| Ñ 0 for the case of a single object B α . This was possible as the current source J 0 is supported away from the object and so H 0 pxq " 1 iωµ 0 ∇ˆE 0 pxq is analytic where the expansion is applied. We extend this to the multiple object case by requiring that J 0 be supported away from B α and introduce the following for n " 1, . . . , N F pnq pxq " 1 2 p∇ zˆE0 pzqqpz pnq qˆpx´z pnq q 1 3 D z p∇ zˆE0 pzqqpz pnq qpx´z pnq qˆpx´z pnq q, ∇ˆF pnq pxq "p∇ zˆE0 pzqqpz pnq q`D z p∇ zˆE0 pzqqpz pnq qpx´z pnq q, so that In other words, ∇ˆF pnq pxq is the first two terms in a Taylor series of iωµ 0 H 0 pxq about z pnq as |x´z pnq | Ñ 0 and so where here and in the following C denotes a generic constant unless otherwise indicated.
Remark 4.1. Higher order Taylor series could be considered (as previously in [13] for the case of a single object) and lead to a more accurate representation of the field in terms of GMPTs. However, in order for such a representation to apply, there will be further implications in the allowable distance between the objects.
The introduction of F pnq pxq motivates the introduction of the following problem: Find w pnq P X α such that where ppB α q pnq q c :" R 3 zpB α q pnq . By the addition of such problems, we havè where w :" ř N n"1 w pnq , w α " w pnq in pB α q pnq and F α " F pnq in pB α q pnq . We also remark that, associated with (13), is the strong form which follows from using Lemma 4.2. For objects pB α q pnq and pB α q pmq with n ‰ m we have that Proof. Introducing ξ pnq " x´z pnq α pnq , which, without loss of generality, we assume the origin to be in B pnq . We set w pnq pxq " α pnq w pnq 0´x´z pnq α pnq¯" α pnq w pnq 0 pξ pnq q and so ∇ xˆw pnq pxq " From the above we have that |w pnq 0 | ď C|ξ pnq |´1}∇ˆE 0 } W 2,8 ppBαq pnq q for sufficiently large |ξ pnq | and so we estimate that |∇ ξˆw pnq 0 | ď C|ξ pnq |´2}∇ˆE 0 } W 2,8 ppBαq pnq q for the same case. Thus, for m ‰ n, Corollary 4.3. Given the description pB α q pnq " α pnq B pnq`zpnq , we are free to configure B pnq in different ways provided that the origin lies at a point in B pnq (similarly with pB α q pmq " α pmq B pmq`zpmq ) . Thus |z pmq´zpnq | will be smallest when the origin lies in the boundaries of the objects, as illustrated in Figure 4. Requiring that |z pmq´zpnq | " min n,m"1,...,N,n‰m |BpB α q pnqB pB α q pmq | ą C ą α max then Lemma 4.2 implies that Figure 4: Illustration to show how each B pnq can be configured differently provided that the origin lies within the object. Consequently d p1q,p2q " |z p1q´zp2q | will be minimum when the objects B p1q and B p2q are configured such that the origin is a suitable point on the boundaries of these objects.
The following Lemma extends Ammari's et al's Lemma 3.2 [2] to the case of N multiple objects, when they are sufficiently well spaced.
Proof. We start by proceeding along the lines presented in [2] and introduce withφ pnq 0 being the solution of an exterior problem in an analogous way toφ 0 in [2]. Using (11) and (14) (and after multiplying by µ 0 ) we can deduce that Also, by application of the Cauchy-Schwartz inequality, we can check that where To bound A 1 and A 2 we have used (12), and applied similar arguments to [2]. The terms A 3 does not appear in the single object case and dictates the minimum spacing for which the bound holds. Requiring that |z pmq´zpnq | " min n,m"1,...,N,n‰m |BB pnq´B B pmq | ą C ą α max and applying Corollary 4.3 then Using (19), (20) and 23) in (18) we find that and, by additionally using }E α´E0´p w pnq`Φpnq q} L 2 pB pnq α q ď α pnq }∇ˆpE α´E0´p w pnqΦ pnq qq} L 2 pB pnq α q ď α max }∇ˆpE α´E0´p w pnq`Φpnq qq} L 2 pB pnq α q , this completes the proof.
By recalling the definition of w Proof. The result immediately follows from Lemma 4.4 and the definition of w pnq 0 .
The expressions for α pnq F pnq pz pnq`αpnq ξ pnq q and w 0 pξ pnq q are obtained by extending in an obvious way the expressions in given in (3.13) and (3.14) in [2] where the latter is now written in terms of pH 0 pz pnq qq i θ pnq i pξ pnq q as well as pD z pH 0 pzqqpz pnq qq ij ψ pnq ij pξ pnq q where θ pnq i pξ pnq q and ψ pnq ij pξ pnq q satisfy the transmission problems and The properties of θ pnq i pξ pnq q and ψ pnq ij pξ pnq q are analogues to the single object case presented in [2].

Integral Representation Formulae
Repeating the proof of Lemma 3.3 in [2] for the multiple object case, it extends in an obvious way to . . Y D pN q be the union of N bounded domains each with Lipschitz boundaries Γ pnq D whose outer normal is n. For any E P H´1p curl; R 3 zDq satisfying ∇ˆ∇ˆE " 0, ∇¨E " 0 in R 3 zD, we have, for any x P R 3 zD In a similar way, repeating the proof of their Lemma 3.4 for multiple objects it extends in an obvious way to

Asymptotic Formulae
Theorem 3.2 in [2] presents the leading order term in asymptotic expansion for pH α´H 0 qpxq for a single inclusion B α as α Ñ 0. In the case of multiple objects that are sufficiently well spaced this extends to Theorem 4.8. For a collection of N objects such that ν pnq is order one, α pnq is small and min n,m"1,...,N,n‰m |BpB α q pnq´B pB α q pmq | ą C ą α max then for x away from B α we have where θ pnq i is the solution of (24) and uniformly in x in any compact set away from B α .
Proof. The proof uses as its starting point Lemma 4.7 and considers each object pB α q pnq in turn. It applies very similar arguments to the proof of Theorem 3.2 in [2] except Theorem 4.5 is used in place of their Theorem 3.1, (22) is used in place of their (3.6) and note that by integration by parts. Furthermore, to recover the negative sign in the first term in (26), we have used as α pnq Ñ 0. Theorem 3.2 in [2] mistakingly uses ∇ x Gpx, αξ`zq " ∇ x Gpx, zq`αD 2 x Gpx, zqξÒ pα 2 q as α Ñ 0, which leads to the wrong sign in their first term, as previously reported for the single homogenous object case in [10].

Results for the Proof of Theorem 3.3
Recall that in this case the object is inhomogeneous and is arranged as B α :" Ť N n"1 B pnq α " α Ť N n"1 B pnq`z " αB`z where α is a single small scaling parameter and z a single translation.

Elimination of the Current Source
The results presented in Section 4.1 hold also in the case when the object is inhomogeneous except the subscript α is replaced by α.

Energy Estimates
For an inhomogeneous object, we proceed along similar lines as [2] and introduce a single vector field F pxq whose curl is such that it is equal to the first two terms of a Taylor series of iωµ 0 H 0 pxq expanded about z as |x´z| Ñ 0 The introduction of F pxq motivates the introduction of the following problem: Find w PX α such that The following Lemma extends Ammari's et al's Lemma 3.2 [2] to the case of an inhomogeneous object.
Lemma 5.1. For an inhomogeneous object B α , there exists a constant C such that for m " 1, . . . , N .
Proof. The result immediately follows from Lemma 5.1 and the definition of w 0 .
The properties of θ i pξq and ψ ij pξq are analogues to the homogeneous object case presented in [2].

Integral Representation Formulae
The integral representation formulae presented in Section 4.3 only require pB α q pnq to be replaced by B pnq α and H α to be replaced by H α for an inhomogeneous object.

Asymptotic Formulae
Theorem 3.2 in [2] presents the leading order term in asymptotic expansion for pH α´H 0 qpxq for a single homogenous inclusion B α " αB`z as α Ñ 0. In the case of an inhomogeneous inclusion this becomes Theorem 5.3. For an inhomogeneous object B α such that ν pnq is order one and α is small then for x away from B α , we have where θ i is the solution of (36) and uniformly in x in any compact set away from B α .
Proof. The proof uses as its starting point Lemma 4.7 and considers each region B pnq α in turn. It applies similar arguments to the proof of their Theorem 3.2 except that our Theorem 5.2 is used in place of their Theorem 3.1 and our (34) instead of their (36). Furthermore, note that by summing contributions, we have that by application of integration by parts and, in a similar manner to the proof of Theorem 4.8, we use to give the correct negative sign in the first term of (38).

Numerical Examples and Algorithms for Object Localisation and Identification
In this section we consider an illustrative numerical application of the asymptotic formulae (5) and (8), numerical examples of the frequency spectra of the MPT coefficients and propose algorithms for multiple object localisation and inhomogeneous object identification as extensions of those in [3].

Numerical Illustration of Asymptotic Formulae for pH α´H 0 qpxq
To illustrate the results in Theorems 3.1 and 3.3, comparisons of pH α´H 0 qpxq 3 will be undertaken with a finite element method (FEM) solver [15] for multiple objects and for inhomogeneous objects. We first show comparisons for two spheres, then comparisons for two tetrahedra followed by comparisons for an inhomogeneous parallelepiped.

Two Spheres
We first consider the situation of two spheres pB α q p1q and pB α q p2q . These objects are defined as pB α q p1q :" which means that the radii of the objects are α p1q and α p2q , respectively. Setting B " B p1q " B p2q to be a sphere of unit radius placed at the origin then are the location of the centroids of the physical objects B p1q α and B p2q α , respectively. Thus, the objects pB α q pnq , n " 1, 2, are centered about the origin with min |BpB α q p1q´B pB α q p2q | " αd. The material properties of the spheres are σ p1q " σ p2q " 5.66ˆ10 7 S/m, µ p1q " µ p2q " µ 0 , we use ω " 133.5rad/s and the object sizes are chosen as α " α p1q " α p2q " 0.01m and hence Mrα p1q B p1q s " Mrα p2q B p2q s, independent of their separation, which will be used in Theorem 3.1. For closely spaced objects we expect Theorem 3.3 to be applicable and in this case we set and z " 0. Note that in this case, M rαBs must be recomputed for each new d.
Comparisons of pH α´H 0 qpxq obtained from the asymptotic formulae (5) and (8) in Theorems 3.1 and 3.3 as well as a full FEM solution are made in Figure 5 for d " 0.2 and d " 2 along three different coordinates axes. To ensure the tensor coefficients were calculated accurately, a p " 3 edge element discretisation and an unstructured mesh of 6 581 tetrahedra is used for computing Mrα p1q B p1q s " Mrα p2q B p2q s and meshes of 8 950 and 11 940 unstructured tetrahedral elements are used for computing M rαBs for d " 0.2 and d " 2, respectively. In addition, curved elements with a quadratic geometry resolution are used for representing the curved surfaces of the spheres. For these, and all subsequent examples, the artificial truncation boundary was set to be 100|B|. To ensure an accurate representation of pH α´H 0 qpxq for the FEM solver, the same discretisation, suitably scaled, as used for M rαBs is employed. For the closely spaced objects, with d " 0.2, we observe good agreement between Theorem 3.3 and the FEM solution in Figure 5, with all three results tending to the same result for sufficiently large |x|. The improvement for larger |x| is expected as the asymptotic formulae (5) and (8) are valid for x away from B α " B α . For objects positioned further apart, with d " 2, we observe that the agreement between Theorem 3.1 and the FEM solution is best. This agrees with what our theory predicts, since, for d " 2, min |BpB α q p1q´B pB α q p2q | " 2α ą α max and so this theorem applies.
Comparisons of pH α´H 0 qpxq for this case are made in Figure 7 for d " 0.2 and d " 2 along three different coordinates axes. To ensure the tensor coefficients are calculated accurately, a p " 3 edge element discretisation and unstructured meshes of 15 617 and 15 488 tetrahedra are used for computing Mrα p1q B p1q s and Mrα p2q B p2q s 4 , respectively, and meshes of 15 837 and 22 045 unstructured tetrahedral elements are used for computing M rαBs for d " 0.2 and d " 2, respectively. To ensure an accurate representation of pH α´H 0 qpxq for the FEM solver, the same discretisation, suitably scaled, as used for M rαBs is employed. As in Section 6.1.1, we observe good agreement between Theorem 3.3 and the FEM solution for the closely spaced objects in Figure 7, with all three results tending to the same result for sufficiently large |x|. For objects positioned further apart, with d " 2, we observe that the agreement between Theorem 3.1 and the FEM solution is again best, which again agrees with what our theory predicts, since, for d " 2, min |BpB α q p1q´B pB α q p2q | " 2α ą α max and so this theorem applies.

Inhomogeneous Parallelepiped
In this section, an inhomogeneous parallelepiped r´1, 0sˆr0, 1sˆr0, 1s Figure 5: Comparison of pH α´H 0 qpxq using the asymptotic expansions (5) and (8) in Theorems 3.1 and 3.1 as well as a FEM solution: along the three coordinate axes for two spheres with different separations αd. Figure 6: Two tetrahedra pB α q p1q and pB α q p2q with min |BpB α q p1q´B pB α q p2q | " αd.
is considered. The material parameters of pB α q p1q and pB α q p2q are µ p1q " µ 0 , σ p1q " 7.371 0 6 S/m, and µ p2q " 5.5µ 0 , σ p1q " 1ˆ10 6 S/m, respectively. Comparisons of pH α´H 0 qpxq obtained from using the asymptotic expansion (8) in Theorem 3.3 and a full FEM solution are made in Figure 8 along three different coordinates axes. To ensure the tensor coefficients are calculated accurately, a p " 3 edge element discretisation and an unstructured mesh of 13 121 tetrahedra are used for computing M rαBs. To ensure an accurate representation of pH α´H 0 qpxq for the FEM solver, the same discretisation, suitably scaled, as used for M rαBs is employed. We observe a good agreement between Theorem 3.3 and the FEM solution for sufficiently large |x|.

Frequency Spectra
The frequency response of the coefficients of MrαBs for a range of single homogeneous objects has been presented in [11,12] where the real part was observed to be sigmoid with respect to log ω and the imaginary part had a distinctive single maxima. Rather than consider the coefficients, it is in fact better to split MrαBs in to the real part RepMrαBsq and an imaginary part ImpMrαBsq, which are both real symmetric rank 2 tensors, and to compute the eigenvalues of these. Indeed, many of the objects previously considered had rotational and/or reflection symmetries such that the eigenvalues coincide with the real and imaginary parts of the diagonal coefficients. A theoretical investigation of RrαBs " RepMrαBs´N 0 rαBsq and IrαBs " ImpMrαBsŃ 0 rαBsq " ImpMrαBsq, where N 0 rαBs corresponds to the real symmetric rank 2 tensor describing the limiting response in the case of ω Ñ 0, and agrees with the Póyla-Szeö tensor for a homogenous permeable object, has been undertaken in [14]. In this we prove results on the eigenvalues of these tensors. Now considering MrαBs for an inhomogeneous object B α , the coefficients of N 0 rαBs are given by N 0 ij rαBs :" where x " x 3 e 3  and we have shown that for 0 ď ω ă 8, the eigenvalues of RrαBs " RepMrαBs´N 0 rαBsq and IrαBs " ImpMrαBs´N 0 rαBsq " ImpMrαBsq have the properties λpRrαBsq ď 0 and λpIrαBsq ě 0 (this also applies to a homogenous objects where B α reduces to B α ) [14].
To illustrate how the behaviour of λpRrαBsq and λpIrαBsq changes for an inhomogeneous object, we consider the geometry of the parallelepiped described in Section 6.1.3 placed at the origin so that B α " B p1q α Y B p2q α " αpB p1q Y B p2q q " αB with α " 0.01m. Note that, although B p1q α and B p2q α have different properties, the object B still reflectional symmetries in the e 1 and e 3 axes and a π{2 rotational symmetry about e 1 so that the independent coefficients of MrαBs are M 11 and M 22 " M 33 (and hence R 11 , R 22 " R 33 are the independent coefficients of RrαBs and I 11 , I 22 " I 33 are the independent coefficients of IrαBs). In Figure 9, we show the computed results for λpRrαBsq and λpIrαBsq for the case where σ p2q " 100σ p1q " 1ˆ10 8 S/m and µ p1q " µ p2q " µ 0 and, in Figure 10, we show the corresponding result for σ p1q " σ p2q " 1ˆ10 6 S/m and µ p2q " 10µ p1q " 10µ 0 . For this we use similar discretisations to those stated previously. In the former case N 0 rαBs vanishes, but not in the latter case. We observe, in Figure 9, that although λ i pRrαBsq, i " 1, . . . , 3 are still monotonically decreasing with log f , it is no longer sigmoid for an inhomogeneous object with varying σ and constant µ and has multiple non-stationary inflection points. Furthermore, rather than a single maxima, λ i pIrαBsq, i " 1, . . . , 3 has two distinct local maxima. However, the results shown in Figure 10 illustrate for an inhomogeneous object with varying µ and constant σ, λ i pRrαBsq, i " 1, . . . , 3, that the behaviour is quite different and, in this case, λ i pRrαBsq, i " 1, . . . , 3 is still sigmoid and the curves for λ i pIrαBsq, i " 1, . . . , 3 still have a single maxima. In the limiting case of ω Ñ 0, λ i pRepMrαBsqq Ñ λ i pRepN 0 rαBsqq, i " 1, . . . , 3 and, for the latter case with a contrast in µ, the behaviour is as shown in Figure 11, which is quite different to a homogenous object of the same size.
Remark 6.1. The results shown in Figures 9 and 12 indicate that the number of points of nonstationary inflection in λ i pRrαBsq and the number of local maxima in λ i pIrαBsq can potentially be used to determine an upper bound on the number of regions with varying σ that make up an inhomogeneous object B α . Note, that contrasts in σ between the regions making up the inhomogeneous object have deliberately chosen as large in these examples and we acknowledge that, for small contrasts, detecting such peaks would be more challenging.

Object Localisation
The approach described by Ammari, Chen, Chen, Volkov and Wang [3] for a single object localisation using multistatic measurements simplifies given our object characterisation in terms of rank 2 MPT for a single homogenous object and also easily extends to inhomogeneous and multiple objects. Following [3], we assume that there are K receivers at locations r pkq , k " 1, . . . , K, λpRepMrαBsqq Figure 14: Frequency dependence of the eigenvalues of RepMrαBsq : inhomogeneous parallelepiped made up of three cubes with σ p1q " σ p2q " σ p1q " 1ˆ10 6 S/m and µ p3q " 10µ p2q " 100µ p3q " 100µ 0 which are associated with small measurement coils with dipole moment q, and L sources at locations s p q , " 1, . . . , L, which are associated with small exciting coils each with dipole moment p. Then, by measuring the field perturbation described by Theorem 3.1 for N target " N objects in the direction q, this gives rise to the k, th entry of the multistatic response matrix as where, for the purpose of the following, we arrange the coefficients of the rank 2 tensors Mrα pnq B pnq s as 3ˆ3 matrices 5 . Assuming that the data is corrupted by measurement noise and is sampled using Hadamard's technique, as in [3], then the MSR matrix can be written in the form whereW " 1 ? 2 pW`iW q and W is a KˆL matrix with independent and identical Gaussian entries with zero mean and unit variance and S noise is a positive constant. In addition, U is a matrix of size Kˆ3N target U "`U p1q¨¨¨U pN q˘, and U pnq is an Kˆ3 matrix U pnq "¨p D 2 Gpr 1 , z pnq qqq 1¨¨¨p D 2 Gpr 1 , z pnq qqq 3 . . . . . . . . .
The matrix M is of size 3N targetˆ3 N target and is block diagonal in the form M " diagpMrα p1q B p1q s,¨¨¨, Mrα pnq B pnq sq, Proceeding in a similar manner to [3], and defining the linear operator L : then, by dropping the higher order term, the MSR matrix can be approximated as The MUSIC algorithm can then be used to localise the location of the multiple arbitrary shaped targets by using the same imaging functional as proposed in [3] I M U pz s q "˜1 where P is the orthogonal projection onto the right null space of LpMq.
Proposition 6.2. Suppose that U M has full rank. Then LpMq has 3N target non-zero singular values. Furthermore, I M U will have N target peaks at the object locations z " z s .
The ability to recover the N target objects will depend on a number of factors : 1. The number and locations of the measurement and excitor pairs. In practice the number of each will be limited to powers of 4 for practical reasons [2].
2. The noise level, which we define as the reciprocal of the signal to noise ratio in terms of the n`3pn´1q th singular value of A 0 (ordered as S 1 pA 0 q ą S 2 pA 0 q . . .
In [2] and [3] the SNR was based instead on the largest singular value S 1 pA 0 q.

The frequency of excitation.
Remark 6.3. From the examination of the frequency dependence of the coefficients of Mrα pnq B pnq s we have seen that the real and imaginary parts for different objects pB α q pnq " α pnq B pnq`zpnq are different. Moreover, in general, their imaginary components exhibit resonance behaviour at different (possibly multiple) frequencies. Consequently, different objects, in general, correspond to different singular values of A 0 . The presence of multiple objects with the same shape and size, but with different locations, will result in multiplicities of the singular values (in the absence of noise). If only a single frequency is considered, and S noise is chosen based on the largest singular value S 1 pA 0 q, then it will generate a W with Gaussian statistics that are associated with only one of the objects possibly present. If the singular values associated with the other objects are much smaller than S 1 pA 0 q it may be difficult to detect the other objects present. In particular, to locate those objects with smaller MPT coefficients (and hence smaller singular values) at that frequency under consideration.
We explore this through the following experiment. We simulate excitations and measurements taken at regular intervals on the plane r´1, 1sˆr´1, 1sˆt0u such that M " K " 256. The dipole moments are chosen as p " q " e 3 so that the plane of all measurement and excitation coils are parallel to this horizontal surface. With these measurements, the location identification of a coin B p1q α of radius 0.01125m and thickness 3.15ˆ10´3m with σ p1q " 15.9ˆ10 6 S/m and µ p1q " µ 0 and a tetrahedron B p2q α with vertices p5.77ˆ10´3, 0, 0qm, p´2.88, 5, 0qˆ10´3m, p´2.88,´5, 0qˆ10´3m and p0, 0,´8.16ˆ10´3qm and material properties σ p2q " 4.5ˆ10 6 S/m and µ p2q " 1.5µ 0 will be considered. The true locations of these objects are assumed to be z p1q " 0.1e 1`0 .1e 2´0 .5e 3 and z p2q "´0.3e 1`0 .3e 2´0 .5e 3 , respectively. To perform the imaging, noise is added to the simulated A 0 to create A and the image functional I M U is evaluated for different positions z s . To do this, we compute P " I M´WS WS where W S are the first 3N singular vectors of A, which are chosen based on the magnitudes of the singular values and thereby allows us to also predict the number of objects N present. We first consider location identification at a frequency f " 1ˆ10 5 Hz, which is close to the resonance peaks for the two objects, and consider the singular values of A 0 and A in Figure 15. At this frequency, S n pA 0 q, n " 1, 2, 3 are associated with the coin and S n pA 0 q, n " 4, 5, 6 with the tetrahedron. Without noise, A " A 0 and the 6 physical singular values are clearly distinguished, but, by considering a noise level of 1% so that S noise " 0.01S 1 pA 0 q, it is no longer possible to distinguish S n pAq, n " 4, 5, 6 from the noisy singular values. On the other hand, by setting S noise " 0.01S 4 pA 0 q, or even S noise " 0.1S 4 pA 0 q, we can distinguish all 6 singular values from the noise. This means that with S noise " 0.01S 1 pA 0 q we expect to only locate the coin, but with S noise " 0.01S 4 pA 0 q, 0.1S 4 pA 0 q we expect to find both objects. This is confirmed in Figure 16 where we plot I M U on the plane´0.5e 3 . We observe that for S noise " 0.01S 1 pA 0 q we can only locate the coin, for S noise " 0.01S 4 pA 0 q we can locate both the coin and the tetrahedron and even by increasing the noise level to 10% and setting S noise " 0.1S 4 pA 0 q both objects can still be identified. On the other hand, choosing the frequency f " 132Hz, such that S n pA 0 q, n " 1, 2, 3 are associated with the tetrahedron and S n pA 0 q, n " 4, 5, 6 with the coin, Figure 17 shows that the phenomena is reversed, and with a 10% noise level and S noise " 0.1S 4 pA 0 q, only the tetrahedron can be identified at this frequency.

Object Identification
A dictionary-based classification technique for individual object identification has been proposed by Ammari et al. [3] and this easily extends to the identification of multiple inhomogeneous objects. We propose a slight variation on that proposed by Ammari et al., which uses the eigenvalues of the real and imaginary parts of the MPT as a classifier as oposed to its singular values at a range of frequencies. The motivation for this is the richness of the frequency spectra of the eigenvalues, as shown in Section 6.2, and that it provides an increased number of values to classify each object. We also propose a strategy in which objects are put in to canonical form before forming the dictionary. The algorithm comprises of two stages as described below.

Off-line Stage
In the off-line stage, given a set of N candidate candidate objects (which can include both homogenous and inhomogeneous objects), we put them in canonical form pB α q piq " α piq B piq`zpiq , i " 1, . . . , N candidate by ensuring that the origin for ξ piq in B piq coincides with the centre of mass of B piq and the object's size α piq is chosen such that |N 0 rα piq B piq s| " 1 6 where N 0 rα piq B piq s " T rα piq B piq s in the case of a homogenous object and corresponds to the Póyla-Szegö tensor as well as being the characterisation for σ˚" 0 for this object 7 . In the case of an object with homogenous materials, the coefficients of Mrα piq B piq s are computed by solving the transmission problem (7) using finite elements and then applying (6) and, in the case of an inhomogeneous object (10) and (9) are used. In each case, the eigenvalues λ R pα piq B piq , ω j q and λ I pα piq B piq , ω j q are obtained for a range of frequencies ω j and D i "tλ R pα piq B piq , ω j q, λ I pα piq B piq , ω j q, j " 1, . . . , N ω u { max k"1,...,Nω p|λ R pα piq B piq , ω k q|, |λ I pα piq B piq , ω k q|q, forms the ith element of the dictionary D " tD 1 , D 2 , . . . , D N candidate u.

On-line Stage
In an extension to [3], the MPT coefficients for each of the targets pT α q piq , i " 1, . . . , N target can be recovered from the same data used to identify the number and locations of objects. Although, to do so, it is important to ensure that the dipole moments of the coils are chosen such that all the 6N target coefficients can be recovered from the measured data [12]. The coefficients are then the solution of the least squares problem pMrpT α q p1q , ω j s,¨¨¨, MrpT α q pN objects q , ω j sq " arg min M }Apω j q´LpMq}, which is repeated for j " 1, . . . , N ω . Then, for each target pT α q piq , we determinê D i " tλ R ppT α q piq , ω j q, λ I ppT α q piq , ω j qu{ max k"1,...,Nω p|λ R ppT α q piq , ω k q|, |λ I ppT α q piq , ω k q|q, and find the closest match toD i within the dictionary D [3]. Notice the target could also be inhomogeneous in which case pT α q piq is replaced by T piq α . 6 For inhomogeneous objects we require |N 0 rα piq B piq s| " 1 and we replace pB α q piq by B piq α " α piq B piq`zpiq , B piq by B piq as well as ensuring the centre of mass coincides with the centre of mass of B piq " Ť N n"1 B pi,nq . 7 If µ˚" µ 0 we choose the object size by requiring the high conductivity limit to have unit determinant.

Numerical Example
As a challenging object identification example, we consider a dictionary consisting of parallelepipeds described in Section 6.2, which consist of either two regions P 1 :" B " B p1q Ť B p2q with B α " αB " αpB p1q Y B p2q q or three regions P 2 :" B " B p1q Y B p2q Y B p3q with B α " αB " αpB p1q Y B p2q Y B p3q q, and vary the material properties according to the descriptions previously described. We also consider the limiting case where the two (three) regions have the same parameters. The dictionary for these objects is generated according to the off-line stage with ω P 2πp2, 300, 4ˆ10 3 , 5ˆ10 4 , 2ˆ10 5 qrad/s, arbitrarily chosen over the frequency spectrum. For the on-line stage take MrT piq α , ω j s, i " 1, 2, j " 1, . . . , N ω " 5 to be given by considering targets T piq α " αRpP i q where R is an arbitrary rotation adding noise. In Figure 18 we illustrate the algorithms ability to differentiate between these similar objects. The red bars indicate the predicated classification, which is correct for the examples presented (it was also found to be correct for the cases of the other parallelepipeds). We can observe that greatest similarity in terms of the classification is between the two homogeneous parallelepipeds and between the two parallelepipeds with contrasting σ and in each case the classification becomes more challenging as the noise level is increased. By increasing the number of frequencies considered so that N ω " 7 with ω P 2πp2, 300, 41 0 3 , 5ˆ10 4 , 2ˆ10 5 , 3ˆ10 6 , 4ˆ10 7 qrad/s we see in Figure 19 that the certainty of the classification is improved for both noise levels.