Electromagnetic Low-Frequency Dipolar Excitation of Two Metal Spheres in a Conductive Medium

This work concerns the low-frequency interaction of a time-harmonic magnetic dipole, arbitrarily orientated in the three-dimensional space, with two perfectly conducting spheres embedded within a homogeneous conductive medium. In such physical applications, where two bodies are placed near one another, the 3D bispherical geometry fits perfectly. Considering two solid impenetrable metallic obstacles, excited by a magnetic dipole, the scattering boundary value problem is attacked via rigorous low-frequency expansions in terms of integral powers ik , where n ≥ 0, k being the complex wave number of the exterior medium, for the incident, scattered, and total nonaxisymmetric electric and magnetic fields. We deal with the static n 0 and the dynamic n 1, 2, 3 terms of the fields, while for n ≥ 4 the contribution has minor significance. The calculation of the exact solutions, satisfying Laplace’s and Poisson’s differential equations, leads to infinite linear systems, solved approximately within any order of accuracy through a cut-off procedure and via numerical implementation. Thus, we obtain the electromagnetic fields in an analytically compact fashion as infinite series expansions of bispherical eigenfunctions. A simulation is developed in order to investigate the effect of the radii ratio, the relative position of the spheres, and the position of the dipole on the real and imaginary parts of the calculated scattered magnetic field.


Introduction
Several practical applications such as Earth's subsurface electromagnetic probing for mineral exploration, geoelectromagnetism, or other physical cases related to the identification of buried metallic or nonmetallic objects of different shapes and sizes stand in the front line of the scientific research.To this end, the primary motivation of the present work came from several geophysical applications, such as mining prospection 1 and detection of cavities 2 or other underground detections for UneXploded Ordinance also referred as UXO 3, 4 and buried objects 5 .Indeed, in those or similar cases we often face the problem of retrieving anomalies of a certain kind, such as metallic ores, from three-dimensional magnetic fields, collected along a borehole when a low-frequency 6 time harmonic source is placed nearby, usually set at the ground surface of the Earth.So, interesting information concerning orientations, sizes, shapes, conductivity, magnetic, and electric parameters of the anomalies can be recovered.
However, this is not an easy task, since solution of this inverse problem often involves a computationally demanding and time-consuming numerical process and it cannot be performed in robust fashion, unless proper models of the field interaction are available.Also, in practice, users want to measure and to identify the body structure at the same time.For this reason, inverse schemes must be fast enough, which can be achieved only if the solution of the direct electromagnetic scattering problem is fast.Therefore, in order to reach this objective, analytical solutions of such physical models are very useful.
There are many computational codes and numerical methods available to that matter.As examples, an inverse scheme is used either to localize a smooth surface of a threedimensional perfectly conducting object using a boundary integral formulation 7 or to find the parameters of an equivalent object 8 using iterative methods.A numerical implementation via integral equations is illustrated in 9 , whilst several complementary investigations are presented in 10 , where the complexity of the applied methods and the difficulty to obtain practical information on buried structures are obvious.Nevertheless, simple methodologies such as to obtain low-frequency expansions by the manipulation of Debye potentials are described in the classical book 11 and remain useful, while a recent and detailed contribution 12 illustrates this type of approach.
Following the assumption of low frequency, there are cases based on the well-known Green's function 13,14 , which assume that the body is a nonmetallic obstacle of finite conductivity.Then, hybrid means from integral formulations of the fields, conveniently approximated, have been developed in order to identify a penetrable sphere illuminated by a dipolar magnetic field within a conductive host medium via a Mie series expansion 15 .In addition, in 16, 17 the problem dealt with is the general case of penetrable obstacles of ellipsoidal shape 18 and it is tackled via expansions of the Green's function or of the electromagnetic fields inside and outside the body, within the framework of the localized nonlinear approximation.
On the other hand, if the body is considered to be metallic, and hence impenetrable, a recent work 19 for perfectly conducting ellipsoids has successfully dealt with the first term the static one of the low-frequency expansion.Moreover, under the basic aim of the present research activities, which is mineral exploration of the Earth by inductive means, unexploded ordinance investigations and exploration of natural structures like water-filled cavities and other possibly conductive materials in subsoil at shallow depths, useful results can be recovered from the already ample library of scattering by simple shapes e.g., spheroids 20 via analytical methods in books 21,22 .At this point, let us mention a successful numerical approach of the characterization of spheroidal metallic objects using electromagnetic induction 23 and an interesting numerical spheroidal-mode approach of unexploded ordinance inversion under time harmonic excitations in the magnetoquasistatic regime 4 .In both applications, complex inversion algorithms are developed.
Despite the fact that, for satisfactory identification of voluminous buried objects, the use of elaborate computer codes is inevitable, there is always need for analytical methods to describe sufficiently the electromagnetic process under consideration.For example, the low-frequency fully 3D electromagnetic fields, scattered from a perfectly conducting sphere in a conductive medium, excited by a magnetic dipole, have been given in a closed analytical form, followed by a numerical demonstration of the results, in 24 .A mathematical formulation of the electromagnetic induction response to spheroidal anomalies embedded in a weakly conducting host medium has been introduced in 25 .Lately, a similar physical model with the same mathematical analysis as in 24 has been introduced in order to extend the geometrical shape of the buried body to the spheroidal 26 and to the ellipsoidal 27 shapes.However, the difficulty of the solution of every physical problem of different geometries from 24 to 26 and, finally, to 27 is strongly increasing as a result of the different kinds of eigenfunctions 14 of the potentials.
Until now, in the above references, only one-body problems are involved.Nevertheless, two-body problems appear in the literature in several physical areas, whereas the most appropriate coordinate system that fits this particular geometry is the bispherical coordinate system 18 , provided that the two obstacles are being considered spherical and set at short distance.Among other publications, low-frequency scattering by two soft spheres is considered in 28 in acoustics, while in 29 the electric potential and field between two different spheres in an infinite medium supporting an external uniform electric field are analytically calculated by the Green's function method.In 30 an analytical solution is developed in order to evaluate the polarizability potential for two spheres, the authors solving the resulting infinite linear systems with cut-off techniques.Although those techniques converge to the final solution only at a certain level of accuracy, they provide us with a good analytical approximation of the fields.
Inductive electromagnetic means which are currently employed in the exploration of the Earth's subsurface and embedded voluminous bodies often call for an intensive use, primary at the modeling stage and later on at the inversion stage, of analytically demanding tools of field calculation.Hence, the already ample library of scattering by simple shapes using strong analytical methods is open to accept new and useful analytical results concerning this matter.This paper describes how to build a versatile set of mathematical tools in order to infer information on two unknown subsurface bodies, by solving the direct scattering problem and calculating the 3D magnetic and electric fields, scattered off when the bodies are illuminated by a magnetic dipole of arbitrary orientation.This primary source is considered to be time harmonic and operated at low frequencies to deeply penetrate in the conductive ground .High-contrast cases, as in our case, for which the ratio between the conductivity of the bodies and the one of the embedding medium is high, are approximated via the assumption of perfectly impenetrable bodies.
To that purpose and according to the well-known low-frequency scattering theory 6 , we expand the electromagnetic fields incident, as well as scattered in positive integral powers of ik , where k is the complex wave number of the exterior medium at the operation frequency.Then, our scattering problem is transformed into an appropriate succession of coupled boundary value scattering problems of the nonaxisymmetric vector fields at each n of the Rayleigh expansion ik n for n ≥ 0. Those problems are formulated according to secondorder Laplace's partial differential equations with proper perfectly reflecting boundary conditions, as a result of the non penetrable character of the boundary of the metallic objects.
Most of the literature e.g., 24 concerns spherical obstacles, a fact that does not hold true for the case of general buried bodies.Although the most general work concerns the ellipsoidal approximation 27 , the simplest and important case is the spheroid 26 .However, here in this paper, we face the problem of two different buried objects placed near each other and, thus, we are obliged to make use of the best fitting system to the particular physics, that is, the bispherical coordinate system 18 .Consequently, we apply the mathematical tools of this nonaxisymmetric geometry to our scattering model, in order to obtain the most important terms of the low-frequency expansions of the electromagnetic fields, those being the static for n 0 and the dynamic n 1, 2, 3 terms.For n ≥ 4 the contribution of the additional terms is of minor significance and they provide no useful information for the fields.
Every problem for n 0, 1, 2, 3 is attacked in the three-dimensional space of the conductive medium, by applying the appropriate boundary conditions on the surface of the spherical bodies, that is, cancellation of the normal component of the total magnetic field and of the two tangential components of the total electric field total fields result from the summation of the incident and scattered fields .After several cumbersome yet rigorous calculations, we analytically calculate the nonaxisymmetric magnetic and electric scattered fields.In particular, for n 0 the Rayleigh term is obtained in terms of an infinite series expansion, while for n 1 there are no incident fields and, thus, no scattered ones.Emphasis is then given on the calculation of the two next nontrivial terms for n 2 and for n 3 of the aforementioned fields, where those are found in closed form from exact solutions of coupled for n 2, to the one at n 0 here we have to solve a Poisson's equation or uncoupled for n 3 Laplace's equations and given in compact fashion as infinite series expansions.
We must mention that the calculation of the exact solution of the zero-order and of the third-order fields leads to infinite linear systems, solved approximately within any order of accuracy through a cut-off procedure.However, this is not the case with the second-order field, where the difficulty of the Poisson's equation to be solved increases the complexity and, hence, a numerical method is employed to solve the corresponding systems, yielded by the boundary conditions.
The rest of the paper itself is organized as follows.In Section 2, the theoretical basis via an analytic mathematical formulation is sketched, while in Section 3 some interesting information concerning the bispherical coordinate system with its limiting cases, as well as the bispherical harmonic eigenfunctions 14 , is provided.In Section 4, the main results of the direct scattering problem at hand are given and discussed, while the difficulties arisen at each step of our calculations are described.Numerical implementation of our results is considered in Section 5, using the analytical formulae derived, and how to get the desired accuracy is discussed.The results include plots that depict the variation of the real and of the imaginary part of the scattered magnetic field as we move upwards or downwards into the subsurface of the Earth.An outline of our work and future steps follow in Section 6.Finally, for completeness, the necessary mathematical material concerning the associated Legendre functions of the first and of the second kind, as well as useful formulae associated with trigonometric and hyperbolic functions, are collected in the appendix, along with the presentation of various vector identities.

Physical and Mathematical Formulation
The particular physics deals with two different solid impenetrable bodies of surfaces S 1 and S 2 under a magnetic dipole excitation.Specifically, we consider two highly conductive spheres buried in a less conductive, homogeneous and isotropic, nonmagnetic medium with conductivity σ and permeability μ, which in fact takes the value of the permeability of free space μ 0 .Provided that the two spheres lie at a reasonable distance from one another, the bispherical geometry introduced next Section 3 , perfectly describes this two-body system.Otherwise, we would have to solve two similar almost identical problems of identifying a metallic sphere 24 .With imaginary unit i, the wave number of the surrounding medium is assumed to be at a low circular frequency ω and with permittivity ε σ/ω, while the harmonic time dependence exp −iωt of all field quantities is implied from now on.The area of electromagnetic action is confined by the external three-dimensional space V R 3 , considered to be smooth and bounded or not bounded, depending on the physical problem.In what follows, every field is written in terms of the position vector r x 1 x 1 x 2 x 2 x 3 x 3 expressed via the Cartesian basis x κ , κ 1, 2, 3 in Cartesian coordinates x 1 , x 2 , x 3 , this dependence omitted for convenience in writing.As for the illumination, it is a known magnetic dipole m with arbitrary location at r 0 and arbitrary orientation: 2.2 If we introduce the H-notation for the magnetic field and the E-notation for the electric field, then the electromagnetic incident fields H i , E i , produced by the magnetic dipole 2.2 , are scattered by the solid spheres, creating the correspondingly scattered fields H s , E s , while the total magnetic and electric fields H t , E t are given by summation of incident and scattered fields, that is, where we have excluded the singular point r 0 from the scattering domain.The boundary value problem is attacked via low-frequency expansions 6 in terms of powers of ik , that is, for the incident i , scattered s , and total t electromagnetic fields.The Maxwell's equations 6 connect magnetic and electric fields and they are written suitably for our purpose as follows: wherein the gradient operator ∇ operates at r.However, whatever the case may be, it could operate at r 0 , consequently for writing convenience we define ∇ ≡ ∇ r and similarly for the Laplace operator Δ ≡ Δ r , otherwise it will be said so.From the low-frequency assumption, the Maxwell's equations 2.5 can be reduced in the form respectively, while all the fields for every r ∈ V R 3 − {r 0 } are divergence-free, as immediate consequence of 2.6 .
Letting R r − r 0 and R |r − r 0 |, the electromagnetic incident fields generated by the magnetic dipole m assume the forms 6

2.8
where the symbol "⊗" denotes juxtaposition.Some extended algebraic calculations on the incident fields 2.8 , based on the Taylor's expansion of the exponential functions e ikR and on definition 2.1 , yield low-frequency relations as powers of ik for the incident fields, where the first four powers of the fields, that is, the static term for n 0 and the dynamic terms for n 1, 2, 3, are sufficient enough.In detail, the primary incident fields enjoy the expressions

2.9
In view of the gradient differential operator ∇ ≡ ∇ r operated at r and the unit dyadic I, we obtain 11 for the magnetic incident fields, while for the electric incident field.In order to derive the second expressions on the right-hand side of the nontrivial incident fields 2.10 -2.13 , which stand for equivalent and easy-to-handle differential forms of those fields, we have used the fact that as well as the result ∇ ⊗ r I, along with identities A.1 and A.2 .It is obvious that the magnetic terms of order n vary like 1/R 3−n , while the electric ones vary like 1/R 4−n as R increases to infinity.We also observe that for the incident magnetic field the dynamic term for n 1 is not present, while for the incident electric field the only term that survives is the dynamic term for n 2, reflecting exactly the same physical and mathematical treatment to the scattered fields.Hence, for the degree of interest here n 0, 2, 3, we obtain the scattered magnetic field and the scattered electric field where the fields H s 0 , H s 2 , H s 3 , and E s 2 are unknown and must be evaluated.Inserting the wave number k of the embedding medium from 2.1 into relations 2.15 and 2.16 , trivial analysis leads to 2.17 respectively.One easily observes that the electric field is purely imaginary, while the magnetic field is a complex number, indicating that the electromagnetic fields at n 2 H s 2 and E s 2 are more than adequate for the full solutions, since the contribution of the H s 3 stands for a very small correction to both the real and the imaginary parts of the scattered magnetic field 2.17 .As far as the zero-order static term H s 0 is concerned, it provides a very good approximation for the real part of the magnetic field 2.17 .Note that there exist no first-order fields as a consequence of lack of the corresponding incident fields resulting from 2.9 .
Recapitulating, we have to calculate the nontrivial scattered fields by solving four mixed Maxwell's type problems for each n 0, 2, 3. Straightforward calculations on Maxwell's equations 2.6 for x s and elaborate use of identity A.7 result in a set of boundary value problems possibly coupled to one another, from the static one at n 0 to dynamic ones at higher values of n up to n 3, given in terms of the harmonic potentials Φ s 0 , Φ s 2 , and Φ s 3 via where the fields H s 0 , H s 2 , E s 2 , and H s 3 are to be calculated in the scattering region 0 as a direct consequence of the type of the incident fields 2.9 and Maxwell's equations 2.6 .The harmonic potentials inside relationships 2.19 -2.21 satisfy the standard Laplace's equations:

2.22
It is worth mentioning that the inhomogeneous vector Laplace equation 2.20 , coupled with the solution of 2.19 , is actually a Poisson's partial differential equation and provided that the zero-order scattered field H s 0 is obtained, then the second-order scattered field H s 2 can be written as the summation of a general vector harmonic function Φ s 2 and of a particular solution 1/2 rΦ s 0 .The form of the particular solution is a straightforward result of the use of identity A.11 with respect to 2.19 .
The set of problems 2.19 -2.21 with 2.22 has to be solved by using the proper perfectly reflecting boundary conditions on the surfaces of the two metal bodies under identification S 1 and S 2 .Those boundary conditions concern the total fields 2.3 at each order n, where in view of the outward unit normal vector n, they demand cancellation of the normal component of the total magnetic field n • H t 0 and of the tangential component of the electric field n × E t 0 .Hence, combining 2.3 and 2.4 , the boundary value problem's conditions yield

2.23
Summarizing now, we are ready to apply the particular bispherical geometry that fits to the two-sphere case, where the boundary value problems to be solved for the magnetic field are the following ones: the static one for n 0, reduced to a potential problem with Neumann boundary condition, the one for n 2, where the problem is far more complicated due to coupling to the static term, and the one for n 3, where we arrive again at a potential problem with Neumann boundary condition.Finally, the scattered electric field for n 2 is given by the curl of the corresponding magnetic field via 2.20 .

Bispherical Geometry and Harmonic Eigenfunctions
In this section, we specialize the previous analysis to the case of two spherical metallic non penetrable obstacles of radii α 1 and α 2 , whose centers are located on the x 3 -axis at a distance d > α 1 α 2 apart.Those spheres are embedded in an otherwise unbounded and continuous medium scattering area , occupying the three-dimensional space V R 3 − {r 0 }.The first task is to adapt exactly one bispherical coordinate system so as to identify the two given spheres with two specified values of the related coordinate variables 18, 28 .Given a fixed positive number c > 0, which we consider to be the semifocal distance of our system, we define the bispherical coordinates ζ, θ, ϕ 18 see Figure 1 , which are connected to the Cartesian ones x 1 , x 2 , x 3 through the equations where ζ ∈ R specifies the nonintersecting spheres for ζ constant, θ ∈ 0, π denotes the intersecting spheres for θ constant, and ϕ ∈ 0, 2π is the azimuthal angle that describes half-planes meridian planes for ϕ constant.Here we must notice that r x 1 , x 2 , x 3 and r 0 x 10 , x 20 , x 30 .If we demand the sphere of radius α 1 in the lower half-space to correspond to the value ζ −ζ 1 on S 1 and the sphere of radius α 2 in the upper half-space to correspond to the value ζ ζ 2 on S 2 with ζ 1 , ζ 2 > 0, then we need to specify ζ 1 , ζ 2 , and c so as to obtain a unique bispherical coordinate system that fits the system of the given twosphere scattering bodies and this is shown in Figure 1.Trivial analysis that appears among several limiting cases in 28 is based on the characteristics of the geometry of the system and provides the unique relations where the radii of the spheres, as well the semifocal distance of the system, have to be much less than the wavelength λ of the incident field, in order to justify the use of lowfrequency approximations.In terms of mathematical analysis, this restriction is secured via the inequality 2c α 1 α 2 λ.Hence, the actual region of propagation of the scattered wave is the exterior domain of the two spheres The spherical radial distance r in the bispherical domain is given by the expression which implies that the far-field region corresponds to a small neighborhood of ζ, θ 0, 0 , whilst one can collapse to some interesting degenerate cases i.e., the plane case for one sphere when ζ 2 ≡ 0 .
The outward unit normal vector n at the surface of each sphere The sign for ζ > 0 changes in order to secure normal differentiation in the outer direction, while, by the definition of a right-handed system ζ, ϕ, θ , the three coordinate vectors are assumed to be 3.9 The gradient ∇ and the Laplacian Δ differential operators also valid for r r 0 are provided by the expressions 3.11 respectively.Obviously, the unit dyadic I can take the bispherical form while the 3D boundary value problem is adjusted to the type of bispherical geometry introduced here, where the x 3 -axis is the axis of symmetry of the scattering region.Before we proceed further, we introduce the interior regular as ζ → −∞ u i mq and the exterior regular as ζ → ∞ u e mq harmonic eigenfunctions of degree 0, 1, 2, . . .and of order m m 0, 1, 2, . . ., in terms of the associated Legendre functions of the first P m and of the second Q m kind 14 see also the appendix via the formulae respectively.Both are regular on the axis of symmetry, while the angular dependence is given by for every ϕ ∈ 0, 2π ,

3.14
with q denoting the even e or the odd o part of the eigenfunctions, the prime denoting derivations with respect to the argument.Then, every harmonic function u in bispherical geometry is written as follows: In addition, as far as the incident fields 2.10 -2.13 with 2.14 are concerned, we utilize the expansion 14 where ε m 1 for m 0 and ε m 2 for m ≥ 1, while r 0 ζ 0 , θ 0 , ϕ 0 .For convenience to our further calculations, we choose to simplify expansion 3.16 and separate the dependence from the singular point r 0 , by introducing, for ≥ 0, m 0, 1, 2, . . ., and q e, o, the function in terms of definition 3.14 and by virtue of the sign function and then

3.20
Concluding, a key formula in works such as ours involving bispherical eigenexpansions is provided by the uniformly convergent expansion 28 In the following section, to solve explicitly the boundary value problems 2.19 -2.23 , we use many useful recurrence relations for the Legendre and the trigonometric functions, as well as certain identities concerning the differential operators, found in appendix.13 for the incident fields, as well as the unit dyadic representation 3.12 , are properly used, in view of the eigenexpansion 3.20 with 3.17 .The position of the magnetic dipole m at r r 0 is also taken into consideration via 2.2 and 3.1 , while harmonic potentials Φ s 0 , Φ s 2 , and Φ s 3 satisfy 2.22 and follow the previous analysis based on relations 3.13 -3.15 .The scattering bispherical domain r ∈ V R 3 − {r 0 } is given by 3.3 , in which the low-frequency magnetic and the electric fields must be built at each n 0, 2, 3. We recall that, in view of 3.4 , there are no electromagnetic fields inside the bodies.

Bispherical Low-Frequency Electromagnetic Fields
Since the region of observation is between the two spheres, we have to employ both interior and exterior harmonic eigenfunctions 3.13 for the potential problems.Under this aspect and with respect to relation 3.10 along with 3.7 -3.9 , we have to perform long and tedious calculations in order to obtain the electromagnetic fields, by making extensive use of the mathematical material of the appendix.In this section, we present the basic and most interesting steps of them.We begin from the easiest case for n 3, continue to the fields for n 0, and conclude with the most cumbersome case for n 2.
The simplest calculations concern the scattered magnetic field H s 3 , due to the fact that the incident field 2.12 for n 3 is a constant vector.Here, we have to solve the potential boundary value problem 2.21 with the Neumann boundary condition 2.23 on both S 1 and S 2 for n 3, which in terms of n given by 3.6 with 3.7 is Then, the interior and the exterior harmonic structure of the potential Φ s 3 , with the help of definitions 3.13 -3.15 , yields for every r ∈ V R 3 − {r 0 }, where a i mq and a e mq for 0, 1, 2, . . ., m 0, 1, 2, . . ., , and q e, o stand for the constant coefficients to be determined.Thus, in terms of the primary field 2.12 , in view of the unit dyadic 3.12 and taking the three projections of the magnetic dipole onto the Cartesian coordinate system from 2.2 , the boundary condition 4.1 , the gradient operator 3.10 and the unit normal vector 3.7 provides us with the following two boundary expressions on   for ≥ 0 are to be evaluated from the two-set boundary conditions 4.6 and 4.7 , applied on However, conditions 4.6 and 4. 7 are not yet ready to profit from the orthogonality of the Legendre functions of the first kind.Indeed, we have to apply the recurrence relation A.20 both on the term cos θP 1 cos θ , ≥ 1 of the first one and the term cos θP cos θ , ≥ 0 of the second one, which appear in the left-hand side of 4.6 and 4.7 , respectively.Then, we rearrange the indexes properly, and by virtue of orthogonality relation A. 19 , we reach systems of two for every value of ζ s ≡ −ζ 1 , ζ 2 linear algebraic equations in the form of sixterm recurrence equations for every κ 1, 2, 3.These equations can be written in a matrix form if we define specific quantities for κ 1, 2, 3.In effect, if as well as the easy-to-handle functions, independent of κ 1, 2, 3, for ≥ 0, 4.10 followed by the simple notations

corresponding to interior constants and
12 corresponding to exterior constants, where j 1, 2 in definitions 4.9 -4.12 stand for the two different spherical surfaces ζ s ≡ −ζ 1 , ζ 2 , respectively.Here, for κ 1, 2 it is valid for ≥ 1, while for κ 3 we have ≥ 0, though, for convenience, we could insert the 0 case into 4.6 without affecting the condition, since in view of A.22 we readily obtain P 1 0 cos θ ≡ 0. Therefore, as long as we force the extra constants and constant coefficients that appear for 0, κ 1, 2 to vanish, we are able to gather the three cases for κ 1, 2, 3 into one general and concrete case.Consequently, the systems of the linear algebraic equations discussed earlier are written as where in view of 4.8 -4.12 , for every 0, 1, 2, . . .and κ 1, 2, 3, the matrix of the coefficients of the unknowns A κ is the six-diagonal symmetric matrix that reads us and is of infinite dimensions.The matrix of constant terms b κ is given by and the matrix of unknown constant coefficients x κ renders where the symbol " " denotes transpose, so as 4.15 and 4.16 be vectors.At this point, it is important to remain compatible with the boundary conditions 4.6 .Hence, especially for the cases at κ 1, 2 the first two lines and the first two columns of the matrix 4.14 must be omitted, in for ≥ 0 can be obtained up to any order of accuracy through cut-off techniques based on the solution of square systems.Once reminding that for n 3 the scattered electric field vanishes, E s 3 0.The efficiency of the cut-off method is demonstrated in Section 5 of this paper, where numerical results produced by the solution of such systems are shown, accompanied by the necessary discussion about the accuracy reached.
A much more complicated analysis, not only in terms of analytical calculations, but also of numerical implementation, is based on the previous steps.Following the same procedure, we are ready to obtain the scattered field when n 0, that is, the static term H s 0 , though not an easy task as immediate consequence of the complexity of the incident field H i 0 , given by 2.10 .This field admits the double action of the gradient operator at position r / r 0 on the quantity 1/R for R |r − r 0 |.Therefore, we are confronted once more with a potential boundary value problem of the form 2.19 and we also apply the Neumann boundary condition 2.23 on both S 1 and S 2 for n 0, whereas for n defined by 3.6 with 3.7 , it is stated by

4.18
Similarly with the previous analysis, the interior and exterior harmonic potential Φ s 0 , with the help of relations 3. for Yet, the expression of the incident field does not appear easily amenable to further processing and an alternative approach has been followed, which is the key to calculation of H s 0 .To do so, we avoid to apply the operator ∇ ≡ ∇ r twice on 1/R, as indicated by relationship 2.10 , but we first evaluate the inner product ζ • ∇ with respect to relations 3.7 and 3.10 to obtain since the dyadic ∇ ⊗ ∇ 1/R is symmetric, whilst the second part of relation 2.14 helps us to avoid the double derivation by changing the argument of derivation over the quantity 1/R.Thus, 4.21 can be rewritten as which, upon introduction of eigenexpansion 3.20 and with the help of the derivative 3.19 , becomes in terms of functions where, in view of 4.24 -4.28 , the matrix of the coefficients of the unknowns A m is valid for ≥ 0, m 0, 1, 2, . . ., and has the six-diagonal symmetric form    and b e mq for 0, 1, 2, . . ., m 0, 1, 2, . . ., , and q e, o up to any order of accuracy through cut-off techniques.Once the solution of the system is obtained, the scattered magnetic field 4.19 , defined in the domain 3.3 , is given by for every r ∈ V R 3 − {r 0 }, while for n 0 there exists no scattered electric field, E s 0 0. Numerical implementation for the static term 4.33 is found in Section 5, where accuracy of the cut-off technique is discussed.
Let us now concentrate upon the potential problem at n 2, where a very difficult and cumbersome manipulation of the boundary value problem 2.20 with 2.23 results into the dynamic scattered fields H s 2 and E s 2 .There exist two reasons that are responsible for this difficulty.The first one is the coupling of the particular model with the zero-order field H s 0 static term , while the second one refers to the extra electric field E s 2 , which enters our problem with the corresponding additional boundary conditions.However, H s 2 as well as E s 2 terms are of major significance, since those fields provide purely imaginary-valued field components within the conductive medium, as seen from 2.17 , 2.18 and contribute to at least most of the imaginaryvalued quadrature part of the magnetic field H s and to the entire imaginary term of the corresponding electric field E s .Indeed, only the real-valued in-phase part of H s is essentially made of the static contribution H s 0 .The mathematical problem to solve is summarized by 2.20 and 2.23 , which in terms of the normal unit vector ζ provided by expression 3.7 becomes

4.35
where the second equality for the scattered electric field in 4.35 comes from immediate application of identity A.4 .Even though the divergence-free character of E s 2 is obvious, this is not the case for the scattered magnetic field H s 2 , where we have which is a consequence of direct application of identity A.3 onto 4.34 using ∇ • r 3. Result 4.36 stands for the extra condition that must be satisfied in addition with the three one scalar and two components of a vector boundary conditions in 4.34 and 4.35 .As readily shown in relation 4.34 , the scattered magnetic field at order n 2 consists of two sets of functions.For the particular solution of the form 2 −1 rΦ s 0 , it is ensured from identity A.11 as well as from the harmonic character of position vector r and potential Φ s 0 that since ∇ ⊗ r I, where coupling with the n 0 solution is exhibited for the nonharmonic part of the field H s 2 .Hence, from 4.18 and 4.19 and in terms of the already calculated constant coefficients b i mq , b e mq for ≥ 0, m 0, 1, 2, . . ., , and q e, o via the system 4.29 -4.32 we write the potential Φ s 0 as for r ∈ V R 3 − {r 0 }.Additionally, the second set of functions of the dynamic field H s 2 is built up from the harmonic character of Φ s 2 for internal and external domains, that is, , given through the scalar ones for ≥ 0, m 0, 1, 2, . . ., , and q e, o.Thus, according to expansions 4.38 and 4.39 and in terms of the bispherical representation of the position vector for every ζ ∈ R and θ ∈ 0, π , the full solution 4.34 for the scattered magnetic field is expressed as while the scattered electric field 4.35 assumes the form , and e i e mq of scalar constant coefficients for ≥ 0, m 0, 1, 2, . . ., , and q e, o, which have been introduced above, have to be calculated in accordance with the incident field data 2.11 and 2.13 , with the particular solution 2 −1 rΦ s 0 and with the boundary conditions 4.34 and 4.35 , provided that the imposed condition of divergence-free scattered magnetic field 4.36 is ensured.
Obviously, the field problem for n 2 has an enormous additional physical difficulty, which is inherited from corresponding mathematical complexity.
Thus, the analytical procedure to follow for the evaluation of the unknown constant coefficients and considered below only emphasizes the main steps, due to the large amount of calculations to perform.We must obtain one double-set relation from boundary where θ ∈ 0, π , ϕ ∈ 0, 2π , and r 0 ζ 0 , θ 0 , ϕ 0 .We must point out that we have omitted many steps, which contain many unnecessary to be shown complicated analytical manipulations, until the form 4.45 is accomplished and we were forced to define certain analytical quantities frequently used to our calculations.The leading functions of the unknown constant coefficients and the constant functions on the right-hand side of 4.45 are defined as follows for κ 1: x 2 e i e mq x 3 for every ≥ 0, m 0, 1, 2, . . ., , and q e, o, after some cumbersome numerical code elaboration.Consequently, we evaluate the scattered magnetic field 4.42 , which satisfies the extra relation 4.36 and the corresponding scattered electric field 4.43 up to a certain level of numerical accuracy, as seen next.

Numerical Results and Discussion
The bispherical coordinate system yields the appropriate environment for solving scattering problems by two spheres.However, this is true only when the low-frequency hypothesis is considered, since we are faced with boundary value problems of Laplace's equation, which accepts R-separation of variables 18 in the bispherical geometry, while Helmholtz's equation does not.For this reason, problems similar to our case adopt this fitting system to obtain analytical results for the corresponding fields.Nevertheless, it not always easy and feasible to pursue fully analytical solutions in closed forms without computational error, or the derived closed formulae do not have the required accuracy, even in the low-frequency realm.Although the electromagnetic fields in the present investigation have been obtained for n 0, 2, 3 there is no first-order field and higher-order terms are not of substantial interest in a closed analytical form of infinite series in terms of the bispherical harmonic eigenexpansions, they are not given in fully compact fashion.Indeed, the constant coefficients are evaluated up to a certain order of accuracy via cut-off techniques for the case of the zeroorder static term and the third-order field or via a direct numerical solution of the boundary condition for the second-order term .Accuracy is controlled by the sum of terms taken in each case, which inherits the semianalytical terminology and character to our method.
For that reason, we adopt the appropriate bispherical geometry to describe two different perfectly conducting σ body → ∞ spheres of radii α 1 , α 2 and we show numerical results about the scattered magnetic field H s ≡ H given by 2.17 , in view of solutions 4.17 , 4.33 , and 4.42 .The centers of the two spheres are on the x 3 -axis at distance d > α 1 α 2 , while embedded in a homogeneous infinite space of magnetic permeability μ 4π • 10 −7 F/m and of electric conductivity σ 2 • 10 −3 S/m.A magnetic dipole as 2.2 , that is, m m 3 x 3 of strength m 3 4π • 10 3 A • m 2 , is located nearby at x 10 , x 20 , x 30 and illuminates the two spheres at the low frequency of 500 Hz, while the magnetic field is measured along a line a borehole at x 1p , x 2p , x 3p .We provide illustrations for the real and the imaginary parts of the approximated scattered magnetic field in units A/m , where the radii and all distances are in meters m .
The case of two spheres of equal radii is considered in Figure 2, where in accordance with the position of the dipole, the results are symmetric with respect to the line x 3 0, as expected.Two similar cases are then examined, with spheres of different radii placed at reasonable distance Figure 3 or located very close, almost touching Figure 4 .In the first set of graphs we observe a small shift of the lines to the x 3 > 0 space, where the large sphere exists, while in the second set there is no significant change, meaning that even for extreme situations we obtain similar results.
The series solution obtained for the evaluation of H converges relatively fast as → ∞, while it must be pointed out that a few terms of the series expansions ∼ 15−25 are  enough for determining the scattered magnetic field within a very good accuracy.Moreover, the number of the adequate terms to obtain a specific accuracy depends upon the distance of the two spheres.As expected, the closer the spheres are, the more terms are needed for achieving the same order of accuracy and vice versa.In addition, for a given accuracy and a given relative distance of the centers of the two spheres, the number of terms needed to converge depends on their ratio of radii, which means an increase of as this ratio increases.Here, let us notice that the above calculations concern the H field.If we would like to check the numerical behavior of the discrete scattered fields for every order separately, we would realize that the terms needed for convergence of the second-order field, where the numerical procedure followed is different, are less numerous than the corresponding necessary number of terms for the cases the zero-order and third-order fields, where the cut-off method is applied.
Finally, let us remark that the effectiveness of the bispherical system in producing analytic or semianalytic results for low-frequency scattering boundary value problems, could serve to obtain more analytical results as a reference tool for more brute-force numerical codes.

Conclusions
The large amount of vector data, the electromagnetic and geometrical complexity of the Earth, the many configurations of sources and receivers, the uncertainty resulting from datasets containing both the contribution of the primary field observed if the ore bodies were, hypothetically, taken out , and the contribution of the secondary field resulting from the interaction of this primary field with the ore bodies explain the continuous interest of elaborating analytical and numerical methods to solve forward and inverse electromagnetic scattering problems.
One of the main goals of the present investigation is to obtain an analytical approximation of the low-frequency behavior of the magnetic and electric fields scattered by two perfectly conducting spheres located close to each other so as not to be considered as two remote objects , illuminated by a time harmonic magnetic dipole.The bispherical geometry matches the particular physical problem of obstacles of smooth surface and arbitrary proportions.As for the assumption of perfectly conducting spherical bodies it is realistic in view of the high conductivity of most mineral ores, their huge conductivity ratio with their surrounding medium, and the low operation frequencies.Hence, we are confronted with a near-field problem, where planar skin depths are significantly larger than the source-body or the body-sensor distances and only diffusion phenomena occur conduction currents are predominant .
The approximation obtained is a low-frequency expansion of the electromagnetic fields up to the third order, since terms of higher order are not of substantial interest and they do not contribute in much sense to the field.The zero-order static term provides a very good approximation concerning the real part of the magnetic field, while the second-order term contributes to the main behavior of its imaginary part.The third-order term stands for a very small correction to both the real and the imaginary parts of the magnetic field, whereas there exist no first-order fields as a consequence of lack of the corresponding incident fields.As for the contribution to the electric field, it comes only from the second-order term, resulting into a pure imaginary form of this field.
Present investigations confirm that simple models as ours appear reliable when used to model the response of two general spheres to a localized vector source in a homogeneous conductive medium both for low-contrast cases and high-contrast cases.Our devised modeling tools are based on a rigorous low-frequency analysis of the electromagnetic fields, where both their real and imaginary parts are of equivalent significance in the development of a reliable model, while further numerical elaboration or implementation in view of a future inversion scheme has been presented.In view of this aspect, mathematical and computational work is currently in progress in several directions, such as different and more complicated geometries or creation of complicated and time-consuming inversion algorithms.

Figure 1 :
Figure 1: The bispherical coordinate system in the three-dimensional space V R 3 , where the various surfaces depicted correspond to different values of ζ ∈ R.

1 , a i e 2 for ≥ 1 and a i e 3
4.15 one have to set D κ 4.6 and must be set to nil.Now we are ready to solve the infinite system 4.13 -4.16 from which the determination of the unknown constant coefficients a i e

. 30 of
infinite dimensions.For ≥ 0, m 0, 1, 2, . . ., and q e, o the matrix of the constant terms b mq is provided by b

cosh ζ cos θ − 1 ζ cosh ζ cos θ − 1 2 ,
straightforward application of identity A.4 , since ∇ × r 0. Obviously, both electromagnetic fields for n 2 are defined in the domain V R 3 − {r 0 } see 3.3 .In this part of the calculations we are forced to work within the purely bispherical geometry.Thus, we represent the fields 4.42 and 4.43 in bispherical coordinates and in order to accomplish it we have to introduce the bispherical image of the constant coefficients h i e mq, which for the prescribedζ ∈ −ζ 1 , ζ 2 , θ ∈ 0, π and ϕ ∈ 0, 2π is h . . .,, and se, o.In order to derive expression 4.44 , we used the unit normal vectors of the system 3.7 -3.9 inserted into the Cartesian form 4.40 .The three sets c condition 4.34 concerning the ζ-component and two double-set relations from boundary condition 4.35 , which correspond to the θ, ϕ-components.The double-set configuration refers to the two particular impenetrable surfaces of the two spherical bodies at S 1 : ζ ζ s ≡ −ζ 1 and S 2 : ζ ζ s ≡ ζ 2 , which offers, finally, six boundary expressions, three for the interior coefficients and three for the exterior coefficients given in 4.40 .In that sense, we apply the boundary conditions 4.34 and 4.35 to the total electromagnetic fields and in view of 4.41 -4.44 , as well as in view of the second-equality form of the incident fields 2.11 and 2.13 , we reach in cumbersome fashion the following six relations for the constant coefficients: q e,o g mq,κ ζ s , θ, ϕ; r 0 for κ 1, 2, 3, ζ s ≡ −ζ 1 , ζ 2 , 4.45

x 3 (Figure 2 :x 3 ( 6 fFigure 3 :
Figure 2: The real and the imaginary parts of the approximated low-frequency scattered magnetic field H for the orders n 0, 2, 3.The two spheres have equal radii α 1 α 2 25 m and their centers are at d 75 m.The magnetic dipole is located at x 10 , x 20 , x 30 200 m, 0 m, 0 m and the measurement line probe is placed at x 1p , x 2p , x 3p 0 m, 200 m, −400, 400 m , obtaining the field on the x 3 -axis.

x 3 (Figure 4 :
Figure 4: The real and the imaginary parts of the approximated low-frequency scattered magnetic field H for the orders n 0, 2, 3.The two spheres have different radii α 1 25 m and α 2 50 m, while their centers are very close d 76 m .The magnetic dipole is located at x 10 , x 20 , x 30 200 m, 0 m, 0 m and the measurement line probe is placed at x 1p , x 2p , x 3p 0 m, 200 m, −400, 400 m , obtaining the field on the x 3 -axis.
2 are given by 2.23 and fit the aforementioned boundary value problems.Since the boundary conditions will be considered upon those two surfaces, it is convenient to introduce the symbol ζ s ≡ −ζ 1 , ζ 2 , which corresponds to ζ s −ζ 1 when we are on S 1 and to ζ s ζ 2 when we are on S 2 .The expressions 2.10 -2.
−1and take orthogonality arguments for the set {cos ϕ, sin ϕ, 1} of trigonometric functions, in order to recover three sets for every κ 1, 2, 3 of two equations for ζ s ≡ −ζ 1 , ζ 2 .Once done, we recall the uniformly convergence expansion 3.21 and we use A.21 so as to obtain o a i mq sinh ζ s 2 1 cosh ζ s − cos θ e 1/2 ζ s a e mq sinh ζ s − 2 1 cosh ζ s − cos θ e − 1/2 ζ s P m cos θ f mq ϕ for ≥ 0 do not appear in expansion 4.2 , since sin 0ϕ 0. In view of 4.4 and definitions 4.5 , we multiply condition 4.3 by 2c cosh ζ s − cos θ 2 |ζ s | P cos θ ζ 2 , 4.8 then in terms of 4.8 on surfaces S 1 at ζ s ≡ −ζ 1 and S 2 at ζ s ≡ ζ 2 , we introduce 19e aid of 3.10 , 3.18 -3.19 taken at r r 0 , while the magnetic dipole m decomposes as shown in 2.2 .Hence, we have achieved the reduction of the difficulty of this boundary condition 4.18 by using this technique.However, if we combine 4.20 and 4.23 , in view of 4.18 , we observe that imposing orthogonality arguments is impossible, since we are obliged to use the recurrence relation A.20 , in order to manipulate properly the terms cos θP m cos θ for ≥ 0 and m 0, 1, 2, ..., , which appear inside both the expressions 4.20 and 4.23 .So, we substitute the modified relationships into condition 4.18 , rearrange the indexes properly, and upon enforcing the orthogonality relation A.19, we similarly to the previous problem for H s 3 reach to systems of two for every value of ζ s ≡ −ζ 1 , ζ 2 linear algebraic six-term recurrence equations.These equations can be expressed in a matrix form if we define specific quantities in order to reduce lengthy relations.Using the sign function 3.18 for the specific cases on surfaces S 1 at ζ s ≡ −ζ 1 and S 2 at ζ s ≡ ζ 2 , while conveniently implying j 1, 2 in the forthcoming definitions in order to identify the two different spherical surfaces ζ s ≡ −ζ 1 , ζ 2 , we introduce s ; r 0 and sign function sgn ζ 0 − ζ s are provided by relation 3.17 and 3.18 , respectively, on surfaces ζ ζ s .The gradient ∇ r 0 ρ mq ζ s ; r 0 is a known quantity, which can be determined with • • • , 4.32 where the symbol " " denotes transpose, so as 4.31 and 4.32 be vectors.So again, we can solve this infinite system 4.29 -4.32 and obtain the unknown constant coefficients b All definitions 4.46 -4.48 hold true for every ζ s ≡ −ζ 1 , ζ 2 with ζ 0 / ζ s , θ ∈ 0, π , and ϕ ∈ 0, 2π , the unit normal vectors ζ ζ s , θ, ϕ , θ ζ s , θ, ϕ , and ϕ ϕ are provided via expressions 3.7 -3.9 , the ad hoc function ρ mq ζ s; r 0 admits the form 3.17 with 3.14 and 3.18 , the magnetic dipole source m is given by 2.2 , and the constant coefficients b i mq and b Even though it would be appropriate to construct this set of boundary conditions in a suitable form by using extensively the pertinent recurrence relations A.20 -A.27 , it seems this effort is futile, since it is impossible by the nature of the mathematics used here to reach a form const. P m cos θ f mq ϕ for 4.45 .That is, we cannot proceed in the same manner as in the previous n 3 and n 0 cases.So, we have handled numerically the boundary conditions 4.45 , in order to obtain the six scalar constant coefficients