## Abstract

The low-frequency electromagnetic scattering problem of a near-surface hollow spherical inclusion is solved by means of a modal approach. The primary field is excited using a current-carrying coil. The presented solution addresses the full coupling between the two interfaces of the geometry in a rigorous way. There is a direct physical interpretation of the different occurring terms in the formal solution, which provides a deeper understanding of the underlying phenomena. The calculation is very fast, which makes the proposed model suitable for use with parametric inversion algorithms.

## 1. Introduction

Low-frequency scattering by spherical inclusions embedded in a conducting medium can be considered as a canonical problem, met in a wide range of applications which range from eddy-current testing (ECT) of conducting materials to light scattering from particles in optics and orebody detection in geophysics. On a second level, the calculation of the scattering response from such canonical inclusions constitutes the first step for addressing the more complex inverse problem of near-surface buried object detection and dimensioning, problem of great importance for geophysical applications.

The scattering from spherical objects above or on planar substrates has been extensively studied in the domain of optics, owing to its theoretical interest and also its particular importance for studying surface phenomena. The first general solution for the electromagnetic scattering of a plane monochromatic wave by a homogeneous sphere in a homogeneous infinite medium was derived in the classical paper of Mie [1]. Using an ingenious combination of Debye and Hertz potentials, Bobbert & Vlieger [2] constructed an analytical solution for the problem of light scattering by a sphere in a non-absorbing medium placed on a plane substrate. The derivation of the formal solution is achieved without approximating assumptions, yet explicit expressions are given for a number of limiting cases such as a perfectly conducting substrate, a small radius as well as a far field region. Videen [3] treats the same problem, but it proceeds to the approximation that the scattered fields emanating from the sphere impinge upon the substrate along the normal direction. Further attempts tackling the same or similar problems are also reported in the literature [4–6].

In all the aforementioned studies, the sphere or the substrate may be of lossy materials, the treatment of highly conducting media is yet not the primary concern of the analysis. Low-frequency scattering by spherical objects consisting of, or embedded in, a conductive medium is a less studied problem (speaking always of analytical or semi-analytical treatment), to the best of the author's knowledge. Vafeas *et al.* [7] have presented an analytical solution for the scattering by a perfectly conducting sphere in a conductive medium using the low-frequency approximation. The eddy-current flow in a conducting sphere induced by an arbitrary current source in the air has been calculated by Theodoulidis *et al.* [8] based on the second-order vector potential formulation (SOVP). Yet, like in the case of optical scattering, the addition of a second boundary of different geometry, such as the planar surface of a half-space, complicates the problem substantially.

In the particular case of low-frequency excitation and the presence of conducting materials, the diffusive nature of the scattering problem allows us to consider the field negligible at a (reasonably long) distance from the excitation, and hence to truncate the solution domain using a simple Dirichlet or Neumann condition, and solve the problem by means of modal methods. This technique, referred to in the literature as truncated region eigenfunction expansion (TREE), has been applied for the solution of a number of important induction problems such as the eddy-current probe interaction with corner discontinuities [9] or boreholes [10,11].

Based on this technique and applying conversion relations between rectangular and cylindrical scalar wave functions, the problem of a cylindrical infinitely long pipe embedded in a conducting half-space is tackled by Skarlatos & Theodoulidis [12]. A key point in this analysis is the exploitation of the common (translational) symmetry axis between the rectangular and the cylindrical coordinate systems used for the treatment of the planar and the cylindrical interface, respectively. In this work, a similar technique is used in order to address the low-scattering problem by a hollow spherical inclusion in a conducting half-space. There is a substantial difference however between this and the previous geometry: here the direction of invariance is the azimuthal angle, thus the use of scalar wave functions is not accessible. This difficulty is recognized by Bobbert & Vlieger [2], who proceed therefore to use a complicated mixing of Debye and Hertz potentials, in order to perform the sought decomposition into scalar problems for the rectangular and spherical coordinate system. Here, in contrast to that work, an analysis similar to Skarlatos & Theodoulidis [12] is carried out, where, however, direct use of vector wave functions in the two coordinate systems is made. The conversion between the wave function of the two systems is based upon the relations presented by Han *et al.* [13]. The presented approach is rigorous, in the sense that no approximation is made apart from that of vanishing fields at far distances from the source.

The paper is organized as follows. After the problem posing and the definition of the vector wave functions in the cylindrical and the spherical coordinate system in §2, the formal solution is derived immediately after (§3). Application of the continuity relations at the interfaces provides the expressions for the development coefficients of the solution in §4. The paper closes with the results section and the conclusions. An alternative expansion, slightly different from that of Han *et al.* [13] for the conversion relations from spherical to cylindrical vector wave functions, is derived in appendix. The derived expansion involves integration over a real-valued wavenumber instead of a real-valued spectral angle as in Han *et al.* [13], which offers increased flexibility when dealing with conductive media.

## 2. Mathematical formulation of the problem

The configuration of the considered problem is depicted in figure 1. A spherical inclusion of radius *R* is embedded in a conducting, non-magnetic half-space with conductivity *σ*. The distance between the centre of the sphere and the half-space interface is *z*_{0}. The material inside the inclusion as well as above the half-space is air.

Eddy-current flow is induced in the half-space by means of an air-cored cylindrical coil, whose axis is normal to the interface. The coil is allowed to move in parallel to the interface at a constant lift-off. Arbitrary current configurations and complex scan trajectories are also allowed by the model, they will not be considered here, however, for the sake of simplicity. The analysis is carried out in the harmonic regime, i.e. the time dependence of the coil current has the form e^{iωt}, *ω* being the radial frequency. The frequency *f*=*ω*/2*π* is considered to be low enough in order for the quasi-static approximation to be valid.

The air half-space above the conductor will be referred to in the text as region 1, the conducting half-space will be named region 2, whereas the index *e* will be reserved for reference to the inclusion volume. We fix the origin of our reference frame to the centre of the spherical inclusion, the *z*-axis being normal to the half-space planar interface as shown in figure 1. Given the distance between the centre of the inclusion and the half-space interface, this means that the latter intersects the *z*-axis at *z*=*z*_{0}. The solution domain is truncated at *ρ*_{L} distance from the *z*-axis, and a perfectly magnetic conductor (PMC) condition is assumed on the resulting cylindrical boundary (that is *B*_{ϕ,z}|_{ρ=ρL}=0).

The source-free induction problem in air is amenable to a magnetostatic formulation via a magnetic scalar potential, i.e. the magnetic flux density can be written as
*Φ*_{1} stands for the potential in air region over the half-space, and *ψ*_{e} is the potential in the inclusion. Both satisfy the Laplace equation
*Φ*_{1} above the half-space comprises the free-space solution *Φ*_{s} ^{1} (i.e. in the absence of the conductor) and the reflection from the half-space *Φ*_{d}

The magnetic induction inside the conducting medium can be written as a sum of two terms, each one associated with one interface. Both are derived by a doublet of scalar potentials [14,15], namely
**c** is a constant pivot vector. The scalar potentials *W*_{a,b} are met in the literature under different names, e.g. Debye, Hertz or second-order potentials; however, all definitions describe more or less the same quantities [15]; in this work, we follow the conventions associated with second-order potential definition. Both satisfy the homogeneous Helmholtz equation
*k*^{2}=i*ωμ*_{0}*σ*≈i*egaμ*_{0}*σ*, *ε*_{0} and *μ*_{0} being the dielectric permittivity and the magnetic permeability of free space, respectively.

Observing that the solution in this part of the structure will have to satisfy the continuity relations on both the planar and the cylindrical surfaces of the problem, we choose to express the total solution as the superposition of two terms, each one being associated with an interface, similarly to Skarlatos & Theodoulidis [12]. Equation (2.4) then can be modified as
**c**:=**e**_{z} and **c**:=**r**, with **e**_{z} being the unit vector along the *z*-direction and **r** the position vector. Recall that we are free to express the solution as any linear combination of partial solutions we wish, under the condition that they will span the corresponding solution space, and they will satisfy the continuity relations across the surfaces of the geometry.

### (a) Vector wave functions

The analysis is facilitated if we make use of the vector wave functions. Setting **c**:=**e**_{z} as a pivot vector, we define the irrotational **L** and the solenoidal vector wave functions **M**,**N** for the cylindrical system as follows
**M** and **N** satisfy the vector wave equation, that is

With **c**:=**r**, **r** being the position vector, we obtain similar definition relations for the spherical wavevector functions **l**,**m** and **n**.

In all the above relations, the following notation convention has been tacitly adopted: all variables associated with the cylindrical coordinate system are denoted using capital letters (i.e. *Φ*_{2},*W*_{a,b},**L**,**M** and **N**), whereas the lower case letters are reserved for variables associated with the spherical coordinate system (namely *ψ*_{e},*w*_{a,b},**l**,**m** and **n**).

## 3. Formal solution

The coordinate system chosen for expressing the solution in region 1 is the cylindrical one, because we need to satisfy the continuity relations on the *z*=*z*_{0} plane. The scalar potentials of (2.3) can then be expanded in Fourier–Bessel series:
*J*_{m}(.) are the cylindrical Bessel functions of the first kind and of order *m*. The *κ*_{n} are determined by the PMC truncation condition on *ρ*_{L}, i.e. *B*_{ϕ,z}|_{ρ=ρL}=0, hence it is

For the solution in the interior of the spherical inclusion, the most appropriate coordinate system to work with is the spherical one. Consequently, defining a spherical coordinate system with origin at the centre of the inclusion, the potential expression in that region can be written as
*m*th order and ℓth degree.

The second-order potentials inside the conducting half-space can be expanded in series of cylindrical and spherical wave functions, as follows
*k*_{n}(.) are the modified spherical Bessel functions of the second kind and of order *n*. Note that only the finite (i.e. vanishing towards *outward*) solutions to the Helmholtz equation have been considered.

Using the vector wave functions definitions and applying the above modal decompositions in (2.1) and (2.6), we obtain for the magnetic induction in the different regions

The vector wave functions in the above relations are obtained from the corresponding scalar ones of the potential expansions after applying the definition relations (2.7)–(2.9) (as well as the respective definitions for the spherical wave functions). Hence, *mn*th mode of the *Φ*_{s} potential, and so on. In addition, note the ‘*o*’ index at the solenoid potential notations, which stands for *outwards* (with respect to the pivot direction **e**_{z} or **r**) in order to distinguish the outwards vanishing solution from the inward one. The necessity of this distinction becomes clear below.

## 4. Matching the modal expressions at the interfaces

### (a) Continuity at the *z*=*z*_{0} plane

In order to address the continuity relations for the magnetic field on the *z*=*z*_{0} plane, we need to express the spherical modes involved in **B**_{2} calculation in terms of cylindrical modes. Following Han *et al.* [13], the expansion of the **m** spherical wave in a basis of cylindrical modes reads
**m**_{mℓ} wave function, whereas the respective relation for the **n**_{mℓ} is easily obtained by applying the properties (2.11), (2.12) (and the respective ones for **m**,**n**), which yields
*v* being defined by *v*^{2}=*κ*^{2}+*k*^{2}. Note that (4.3), (4.4) are slightly different than those in Han *et al.* [13]. By virtue of the truncation relation introduced in (3.3), the above integrals are converted into series, i.e.
*a*_{mℓn},*b*_{mℓn} are related by the continuous ones via the relation
*b*_{mℓn}; *ρ*_{L} is the truncation boundary.

Substituting (4.5), (4.6) into (3.8) yields

From the *H*_{ρ} continuity at *z*=*z*_{0} and taking into account the orthogonality of the angular modes, we obtain for the mode *m* (after term-rearrangement):

In the same fashion, the continuity relations for *H*_{ϕ} and *B*_{z} yield
*J*_{m} in (4.9) complicates the analysis, therefore we use the Bessel function identities

Weighting (4.13) with the *J*_{m±1} over the interval [0,*ρ*_{L}] and taking into account the relation [16]
*δ*_{nn′} is the Kronecker delta, we arrive at

In the same fashion, the application of the orthogonality relation to (4.11) yields

### (b) Continuity on the spherical surface

Because the observation is made on the spherical surface, we need to proceed to the opposite transformation than before, namely the cylindrical functions have to be expressed in terms of spherical ones. These conversion relations are given in Han *et al.* [13], and read
**M** wave function, and
**N** one. The expansion coefficients *A*_{mℓ} and *B*_{mℓ} are given by the relations
*α* has been defined in the previous subsection.

Substituting, as above, the wave function expressions (4.17)–(4.18) into (3.8) and applying the continuity relations at the surface of the sphere *r*=*R* for the magnetic field components, we obtain the following relations for the mode *m*; for *B*_{r}:
*H*_{ϕ}:
*H*_{θ}:

We introduce now the following orthogonality relations for the associated Legendre polynomials [17]:
*δ*_{ℓℓ′} is the Kronecker delta, as usual.

Applying (4.24) to (4.21), it yields
*π*] and taking (4.25) into account yields

### (c) Construction of the discrete system

Summarizing the previously derived coefficient relations, we have

Following the usual procedure, we truncate the infinite series in order to produce a finite system of linear equations, whose solution will yield the unknown development coefficients. General rules about the series truncation criteria have been given in previous works [9,10]. Hence, we are led to the following finite algebraic system:
**d**_{d},**c**_{a},**c**_{b},**d**_{a},**d**_{b} and **d**_{e} stand for the coefficient vectors.

Note here that the different scattering terms related to the two interfaces as well as their coupling through the

Finally, the coil impedance variation owing to the presence of the half-space with the inclusion can be calculated in the standard way using the reciprocity theorem, which leads to the well-known expression [15]:

## 5. Results

The presented model has been applied to the inspection situation depicted in figure 1, where spherical inclusions with different radii and a highly conducting half-space are considered. The half-space conductivity is *σ*=35.4 MS m^{−1}. The coil inner and outer radius is 2 and 4 mm, respectively, its length is 1 mm and it is wound with 200 turns. The coil is moved at a constant lift-off equal to 0.2 mm from the half-space interface.

The results are compared against numerical simulations obtained using the three-dimensional finite-element method (FEM) package Comsol 3.5a in figure 2 for two frequencies: 1 and 5 kHz. The first inclusion has 5 mm radius and its ligament from the half-space free surface is 0.5 mm. The second one is smaller with 1 mm radius and a ligament of 0.1 mm.

As a figure of merit, the truncation zone for the TREE solution can be taken equal to about 10 times the radial extent of the coil (provided that the inclusion volume is entirely contained in the domain defined by this limit). Truncation radii greater than the previous one have an almost negligible impact to the accuracy of the results. From the computational point of view, larger *ρ*_{L} values impose a larger number of radial modes to be taken into account, which leads to a cubic increase (due to the full system matrix) of the computational time.

The calculation time for the calculation of a complete scan comprising 60 points and for the above-mentioned truncation limit was at most 5 s for both inclusions and both frequencies. The respective FEM calculation time instead is estimated to 1.5 min per scan point, both times being measured using a standard workstation. It should be pointed out that the computational effort for the semi-analytical solution rises slowly with the number of scan points, because the system matrix is independent of the coil position, the latter affecting only the right-hand-side vector of the linear system. Therefore, using LU-variants for the system inversion, the additional computational cost is merely due to the backwards substitution calculations. This is a well-reported advantage of semi-analytical solutions [9,10].

The reason for the very small calculation times lies also with the rapid convergence of the solution series, which allowed us to consider few modes. More precisely, the number of cylindrical modes did not exceed 80 for *κ* and 25 for *m* for all cases, whereas the maximum number of spherical modes was at most 15. This is a non-surprising remark, because the cylindrical and spherical wave functions bases used for the solution expansion are partial solutions of the given geometry, and thus they are already very close to the final solution.

## 6. Conclusion

A semi-analytical solution to the low-frequency scattering problem for a hollow inclusion in a conducting half-space is derived. The full coupling between the two interfaces is taken into account without using approximations, apart from the one of a vanishing field at sufficiently long distance from the source. The error introduced by this assumption, however, can be made arbitrarily small just by fixing the truncation distance.

The calculation is very fast, which makes the presented model a suitable candidate for use with inversion algorithms.

## Acknowledgements

This work was conducted in the context of the CIVAMONT collaboration (http://www-civa.cea.fr).

## Appendix A. Derivation of the conversion relations from spherical to cylindrical vector wave function

The starting point for the following derivation is the conversion relation from spherical to cylindrical scalar wave functions given in Chang & Mei [18]:
*u*]≥0. Using the (more common for the eddy-current induction problems) *k*^{2} sign convention of (2.5) as well as the identities

We seek a conversion relation of the form^{2}
*a*_{nm} and *b*_{nm} are to be determined.

We form the inner product of (A 8) with **M**_{−m}, and we obtain

We consider each integral separately. We begin with the first integral on the right-hand side. Using the cylindrical vector wave function definition (2.8), we have
*I*_{1})
*V* encloses the positive or negative half-space, i.e. it consists of the *z*=0 interface and an hemisphere with radius ^{3}

In the same way, we obtain for the second integral of the right-hand side of (A 9) (denoted as *I*_{2})
*I*_{1}, we state that the boundary integral vanishes, and for the second, we obtain (taking (2.12) into account)
**M** and **N** are orthogonal to each other.

Let us now consider the left-hand side integral of (A 9) (let us denote it *I*_{3}). Following the same procedure, we have
*c*_{nm} from (A 7) yields

From (A 9), (A 15), (A 17) and (A 25), we obtain our final result for the coefficient *a*_{mn}:

For the calculation of the second coefficient, *b*_{mn}, we return to the general expression for the **m**_{mn} development (A 8), and we apply the curl operator on both sides, which using the properties of the vector wave functions becomes

Projecting (A 27) upon **M**_{−m} as before, we obtain for the left-hand side

Again, the boundary integral is zero, as it can be easily verified, and for the remaining term, we obtain after some manipulations the result
*a*_{mn} expression), we obtain for *b*_{mn}
*c*_{mn}, we arrive at the final result

## Footnotes

↵1 Strictly speaking,

*Φ*_{s}satisfies the homogeneous Laplace equation. Nevertheless, because we are merely interested in the solution just above the half-space interface, we restrict the analysis to the source-free region between the coil and the interface, i.e.*z*_{0}<*z*<*z*_{b},*z*_{b}being the coil base.↵2 Here, the explicit distinction between inwards and outwards evanescent functions is abandoned, because there is no risk of confusion.

↵3 Note that the

*z*-signs in the corresponding integral are combined in the same way thus delivering the same expression, regardless the half-space considered.

- Received April 3, 2014.
- Accepted July 23, 2014.

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