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

Authors: Panayiotis Vafeas; Polycarpos K. Papadopoulos; Dominique Lesselier
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

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)n, where n≥0, k being the complex wave number of the exterior medium, for the incident, scattered, and total non-axisymmetric 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.


---

## Body

## 1. 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 three-dimensional 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 second-order 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 forn=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 Section2, 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.


## 2. Physical and Mathematical Formulation

The particular physics deals with two different solid impenetrable bodies of surfacesS1 and S2 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(2.1)k=iωμσ=ωμσ2(1+i),
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(ℝ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=x1x̂1+x2x̂2+x3x̂3 expressed via the Cartesian basis x̂κ, κ=1,2,3 in Cartesian coordinates (x1,x2,x3), this dependence omitted for convenience in writing. As for the illumination, it is a known magnetic dipole m with arbitrary location at r0 and arbitrary orientation:(2.2)m=∑i=13mix̂i.If we introduce theH-notation for the magnetic field and the E-notation for the electric field, then the electromagnetic incident fields Hi,Ei, produced by the magnetic dipole (2.2), are scattered by the solid spheres, creating the correspondingly scattered fields Hs,Es, while the total magnetic and electric fields Ht,Et are given by summation of incident and scattered fields, that is,(2.3)Ht=Hi+Hs,Et=Ei+Esforr∈V(R3)-{r0},
where we have excluded the singular point r0 from the scattering domain. The boundary value problem is attacked via low-frequency expansions [6] in terms of powers of (ik), that is,(2.4)Hx=∑n=0∞Hnx(ik)n,Ex=∑n=0∞Enx(ik)nforr∈V(R3)-{r0},x=i,s,t,
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:(2.5)∇×Ex=iωμHx,∇×Hx=(-iωε+σ)Ex≅ε≪σ/ωσExwith∇⋅Hx=∇⋅Ex=0,x=i,s,t,
wherein the gradient operator ∇ operates at r. However, whatever the case may be, it could operate at r0, 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(2.6)σ∇×Enx=-Hn-2xforn≥2,∇×Hnx=σEnxforn≥0,x=i,s,t,
respectively, while all the fields for every r∈V(ℝ3)-{r0} are divergence-free,(2.7)∇⋅Hnx=∇⋅Enx=0forn≥0,x=i,s,t,
as immediate consequence of (2.6).LettingR=r-r0 and R=|r-r0|, the electromagnetic incident fields generated by the magnetic dipole m assume the forms [6](2.8)Hi=14π[(k2+ikR-1R2)m-(k2+3ikR-3R2)R⊗R⋅mR2]eikRR,Ei=[ωμk4π(1+ikR)m×RR]eikRR,
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 eikR 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)Hi=[H0i+H2i(ik)2+H3i(ik)3]+◯((ik)4),Ei=[E2i(ik)2]+◯((ik)4).
In view of the gradient differential operator ∇≡∇r (operated at r) and the unit dyadic Ĩ, we obtain (2.10)H0i=m4π⋅(3R⊗RR2-Ĩ)1R3=m4π⋅(∇⊗∇1R),(2.11)H2i=-m4π⋅(R⊗RR2+Ĩ)12R=m4π⋅12(∇1R⊗R-ĨR),(2.12)H3i=m4π⋅(-23Ĩ)
for the magnetic incident fields, while(2.13)E2i=-m4πσ×RR3=m4πσ×∇1R
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(2.14)∇r1R≡∇1R=-RR3,∇r1R≡∇1R=-∇r01R,
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/R3-n, while the electric ones vary like 1/R4-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 heren=0,2,3, we obtain the scattered magnetic field(2.15)Hs=[H0s+H2s(ik)2+H3s(ik)3]+◯((ik)4)forr∈V(R3)-{r0}
and the scattered electric field(2.16)Es=[E2s(ik)2]+◯((ik)4)forr∈V(R3)-{r0},
where the fields H0s,H2s,H3s, and E2s 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)Hs=[H0s+(ωμσ)ωμσ2H3s]+(ωμσ)[ωμσ2H3s-H2s]i+◯((ik)4)forr∈V(R3)-{r0},(2.18)Es=[-(ωμσ)E2s]i+◯((ik)4)forr∈V(R3)-{r0},
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 (H2s and E2s) are more than adequate for the full solutions, since the contribution of the H3s 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 H0s 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 eachn=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 Φ0s,Φ2s, and Φ3s via(2.19)ΔH0s=0⟹H0s=∇Φ0s,since∇⋅H0s=0,∇×H0s=0,(2.20)ΔH2s=H0s⟹H2s=Φ2s+12(rΦ0s),σE2s=∇×H2s,since∇⋅H2s=∇⋅E2s=0,(2.21)ΔH3s=0⟹H3s=∇Φ3s,since∇⋅H3s=0,∇×H3s=0,
where the fields H0s,H2s,E2s, and H3s are to be calculated in the scattering region V(ℝ3)-{r0}, while H1s=E0s=E1s=E3s=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)ΔΦ0s=ΔΦ3s=0,ΔΦ2s=0.
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 H0s is obtained, then the second-order scattered field H2s can be written as the summation of a general vector harmonic function Φ2s and of a particular solution (1/2)(rΦ0s). 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 S1 and S2. 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̂·Ht=0) and of the tangential component of the electric field (n̂×Et=0). Hence, combining (2.3) and (2.4), the boundary value problem’s conditions yield(2.23)n̂⋅(Hni+Hns)=0forn=0,2,3,n̂×(E2i+E2s)=0foreveryr∈S1∪S2.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 forn=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).


## 3. 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 x3-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(ℝ3)-{r0}. 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 (x1,x2,x3) through the equations(3.1)x1=csin⁡⁡θcos⁡⁡φcosh⁡⁡ζ-cos⁡⁡θ,x2=csin⁡⁡θsin⁡⁡φcosh⁡⁡ζ-cos⁡⁡θ,x3=csinh⁡⁡ζcosh⁡⁡ζ-cos⁡⁡θ,
where ζ∈ℝ 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=(x1,x2,x3) and r0=(x10,x20,x30). If we demand the sphere of radius α1 (in the lower half-space) to correspond to the value ζ=-ζ1 on S1 and the sphere of radius α2 (in the upper half-space) to correspond to the value ζ=ζ2 on S2 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 two-sphere 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(3.2)c=d2-(α1-α2)2d2-(α1+α2)22d,ζκ=sinh⁡-1⁡(cακ)forκ=1,2,
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 low-frequency approximations. In terms of mathematical analysis, this restriction is secured via the inequality 2c+α1+α2≪λ.Figure 1
The bispherical coordinate system in the three-dimensional spaceV(ℝ3), where the various surfaces depicted correspond to different values of ζ∈ℝ.Hence, the actual region of propagation of the scattered wave is the exterior domain of the two spheres(3.3)V(R3)-{r0}={(ζ,θ,φ):ζ∈(-ζ1,ζ2),θ∈[0,π],φ∈[0,2π)}
and forms a bispherical shell. Since they are impenetrable, there is no wave propagation inside them, hence, there are no electromagnetic fields produced either in the interior region ζ∈(-∞,-ζ1) of the first sphere S1 (at ζ=-ζ1) or in the interior region ζ∈(ζ2,+∞) of the second sphere S2 (at ζ=ζ2). Therein,(3.4)Hx=Ex=0,x=i,s,tforζ∈(-∞,-ζ1)∪(ζ2,+∞),θ∈[0,π],φ∈[0,2π).
The spherical radial distance r in the bispherical domain is given by the expression(3.5)r=ccosh⁡⁡ζ+cos⁡⁡θcosh⁡⁡ζ-cos⁡⁡θ,wherec>0,ζ∈(-ζ1,ζ2),θ∈[0,π],
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 vectorn̂ at the surface of each sphere ζ=-ζ1 or ζ=ζ2 is provided in the Cartesian basis x̂κ,κ=1,2,3, by the formula(3.6)n̂=|sinh⁡⁡ζ|sin⁡⁡θcos⁡⁡φcosh⁡⁡ζ-cos⁡⁡θx̂1+|sinh⁡⁡ζ|sin⁡⁡θsin⁡⁡φcosh⁡⁡ζ-cos⁡⁡θx̂2+|sinh⁡⁡ζ|(cosh⁡⁡ζcos⁡⁡θ-1)sinh⁡⁡ζ(cosh⁡⁡ζ-cos⁡⁡θ)x̂3={ζ̂,ζ<0,x̂3,ζ=0,-ζ̂,ζ>0.
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.7)ζ̂=-1cosh⁡⁡ζ-cos⁡⁡θ[sinh⁡⁡ζsin⁡⁡θ(cos⁡⁡φx̂1+sin⁡⁡φx̂2)+(cosh⁡⁡ζcos⁡⁡θ-1)x̂3],(3.8)θ̂=1cosh⁡⁡ζ-cos⁡⁡θ[(cosh⁡⁡ζcos⁡⁡θ-1)(cos⁡⁡φx̂1+sin⁡⁡φx̂2)-sinh⁡⁡ζsin⁡⁡θx̂3],(3.9)φ̂=-sin⁡⁡φx̂1+cos⁡⁡φx̂2.The gradient∇ and the Laplacian Δ differential operators (also valid for r=r0) are provided by the expressions(3.10)∇r≡∇=cosh⁡⁡ζ-cos⁡⁡θc[ζ̂∂∂ζ+θ̂∂∂θ+φ̂sin⁡⁡θ∂∂φ],(3.11)Δr≡Δ=(cosh⁡⁡ζ-cos⁡⁡θ)2c2sin⁡⁡θ×{(cosh⁡⁡ζ-cos⁡⁡θ)[∂∂ζ(sin⁡⁡θcosh⁡⁡ζ-cos⁡⁡θ∂∂ζ)+∂∂θ(sin⁡⁡θcosh⁡⁡ζ-cos⁡⁡θ∂∂θ)]+1sin⁡⁡θ∂2∂φ2},
respectively. Obviously, the unit dyadic Ĩ can take the bispherical form(3.12)Ĩ=∑i=13x̂i⊗x̂i=ζ̂⊗ζ̂+θ̂⊗θ̂+φ̂⊗φ̂,
while the 3D boundary value problem is adjusted to the type of bispherical geometry introduced here, where the x3-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(3.13)ul(i)mq=e(l+1/2)ζPlm(cos⁡⁡θ)fmq(φ),ul(e)mq=e-(l+1/2)ζPlm(cos⁡⁡θ)fmq(φ),
respectively. Both are regular on the axis of symmetry, while the angular dependence is given by(3.14)fmq(φ)={cos⁡⁡mφ,q=esin⁡⁡mφ,q=o⟹fmq′(φ)={-msin⁡⁡mφ,q=emcos⁡⁡mφ,q=oforeveryφ∈[0,2π),
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:(3.15)Δu=0⟹u=cosh⁡⁡ζ-cos⁡⁡θ∑l=0∞∑m=0l∑q=e,o[Al(i)mqul(i)mq+Al(e)mqul(e)mq],
where Aℓ(i)mq and Aℓ(e)mq are constant coefficients.In addition, as far as the incident fields (2.10)–(2.13) with (2.14) are concerned, we utilize the expansion [14](3.16)1R≡1|r-r0|=1c(cosh⁡⁡ζ-cos⁡⁡θ)(cosh⁡⁡ζ0-cos⁡⁡θ0)×∑l=0∞∑m=0lεm(l-m)!(l+m)!e-(l+1/2)|ζ-ζ0|Plm(cos⁡⁡θ)Plm(cos⁡⁡θ0)cos⁡⁡[m(φ-φ0)],
where εm=1 for m=0 and εm=2 for m≥1, while r0=(ζ0,θ0,φ0). For convenience to our further calculations, we choose to simplify expansion (3.16) and separate the dependence from the singular point r0, by introducing, for ℓ≥0,m=0,1,2,…,ℓ and q=e,o, the function(3.17)ρlmq(ζ;r0)=εmc(l-m)!(l+m)!cosh⁡⁡ζ0-cos⁡⁡θ0e-[sgn⁡(ζ0-ζ)](l+1/2)ζ0Plm(cos⁡⁡θ0)fmq(φ0),
in terms of definition (3.14) and by virtue of the sign function(3.18)sgn⁡(ζ0-ζ)={-1,ζ0<ζ,0,ζ0=ζ,1,ζ0>ζ,andthensgn⁡(ζ0-ζ)=ζ0-ζ|ζ0-ζ|={-1,ζ=ζ2,1,ζ=-ζ1,sinceζ0≠-ζ1,ζ2,
accompanied by its ζ-derivative(3.19)∂∂ζ[sgn⁡(ζ0-ζ)]=-2δ(ζ0-ζ)={0,ζ0≠ζ,-2,ζ0=ζ,
where δ(ζ0-ζ) is the delta function. Consequently, since |ζ-ζ0|=-[sgn(ζ0-ζ)](ζ-ζ0) and by appropriate combination of (A.24) with (A.25), expansion (3.16) is transformed to(3.20)1R≡1|r-r0|=cosh⁡⁡ζ-cos⁡⁡θ∑l=0∞∑m=0l∑q=e,o[ρlmq(ζ;r0)e[sgn⁡(ζ0-ζ)](l+1/2)ζPlm(cos⁡⁡θ)fmq(φ)].Concluding, a key formula in works such as ours involving bispherical eigenexpansions is provided by the uniformly convergent expansion [28](3.21)1cosh⁡⁡ζ-cos⁡⁡θ=2∑l=0∞e-(l+(1/2))|ζ|Pl(cos⁡⁡θ).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.


## 4. Bispherical Low-Frequency Electromagnetic Fields

Our purpose is to derive closed analytical expressions for the nonvanishing scattered magnetic fieldsH0s,H2s, and H3s, as well as for the scattered electric field E2s, since H1s=E0s=E1s=E3s=0. To achieve it, we must independently solve the problems (2.19) and (2.21) to get H0s and H3s, respectively, and then proceed to the problem (2.20) to get H2s and, thus, E2s, which is much more complicated due to coupling with (2.19). The properly adjusted boundary conditions for the total electromagnetic fields (2.3) on the two surfaces ζ=-ζ1 and ζ=ζ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 S1 and to ζs=ζ2 when we are on S2. The expressions (2.10)–(2.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=r0 is also taken into consideration via (2.2) and (3.1), while harmonic potentials Φ0s,Φ2s, and Φ3s satisfy (2.22) and follow the previous analysis based on relations (3.13)–(3.15). The scattering bispherical domain r∈V(ℝ3)-{r0} 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.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 fieldH3s, 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 S1 and S2 for n=3, which in terms of n̂ given by (3.6) with (3.7) is(4.1)H3s=∇Φ3swithΔΦ3s=0,whereζ̂⋅(H3i+H3s)=0on⁡ζ=ζs≡-ζ1,ζ2.
Then, the interior and the exterior harmonic structure of the potential Φ3s, with the help of definitions (3.13)–(3.15), yields(4.2)H3s=∑l=0∞∑m=0l∑q=e,o∇[cosh⁡⁡ζ-cos⁡⁡θ(al(i)mqe(l+1/2)ζ+al(e)mqe-(l+1/2)ζ)Plm(cos⁡⁡θ)fmq(φ)]
for every r∈V(ℝ3)-{r0}, 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 ζ=ζs≡-ζ1,ζ2:(4.3){cosh⁡⁡ζs-cos⁡⁡θ2c×∑l=0∞∑m=0l∑q=e,o[al(i)mq(sinh⁡⁡ζs+(2l+1)(cosh⁡⁡ζs-cos⁡⁡θ))e(l+1/2)ζs+al(e)mq(sinh⁡⁡ζs-(2l+1)(cosh⁡⁡ζs-cos⁡⁡θ))e-(l+1/2)ζs]Plm(cos⁡⁡θ)fmq(φ)+2314π1cosh⁡⁡ζs-cos⁡⁡θ×[sinh⁡⁡ζssin⁡⁡θ(m1cos⁡⁡φ+m2sin⁡⁡φ)+m3(cosh⁡⁡ζscos⁡⁡θ-1)]∑m=0l}=0.
An immediate consequence of the last summand of relation (4.3), corresponding to incident field, is that(4.4)al(i)mq=al(e)mq=0form=2,3,…,l,l≥2,q=e,o,
while, for convenience, we adopt the following notation for the nonvanishing coefficients:(4.5)al(ie)1e≡al(ie)1,al(ie)1o≡al(ie)2forl≥1,al(ie)0e≡al(ie)3forl≥0.
Note that the coefficients aℓ(ie)0o for ℓ≥0 do not appear in expansion (4.2), since sin0φ=0. In view of (4.4) and definitions (4.5), we multiply condition (4.3) by 2c(coshζs-cosθ)-1 and 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(4.6)∑l=1∞[al(e)κ((2l+1)(cosh⁡⁡ζs-cos⁡⁡θ)-sinh⁡⁡ζs)e-(l+1/2)ζs-al(i)κ((2l+1)(cosh⁡⁡ζs-cos⁡⁡θ)+sinh⁡⁡ζs)e(l+1/2)ζs]Pl1(cos⁡⁡θ)=cmκsinh⁡⁡ζs3πsin⁡⁡θ(cosh⁡⁡ζs-cos⁡⁡θ)3/2=-2cmκsinh⁡⁡ζs3π∂∂θ1cosh⁡⁡ζs-cos⁡⁡θ=-22cmκsinh⁡⁡ζs3π∑l=0∞e-(l+1/2)|ζs|∂∂θPl(cos⁡⁡θ)=22cmκsinh⁡⁡ζs3π∑l=0∞e-(l+1/2)|ζs|Pl1(cos⁡⁡θ)
for κ=1,2. On the other hand, we apply one more time the very useful expansion (3.21), which in combination with equality |ζs|=(sgnζs)ζs and the properly adjusted derivative (3.19) exhibits(4.7)∑l=0∞[al(e)3((2l+1)(cosh⁡⁡ζs-cos⁡⁡θ)-sinh⁡⁡ζs)e-(l+1/2)ζs-al(i)3((2l+1)(cosh⁡⁡ζs-cos⁡⁡θ)+sinh⁡⁡ζs)e(l+1/2)ζs]Pl(cos⁡⁡θ)=cm33π(cosh⁡⁡ζscos⁡⁡θ-1)(cosh⁡⁡ζs-cos⁡⁡θ)3/2=cm33π(sinh⁡2ζs(cosh⁡⁡ζs-cos⁡⁡θ)3/2-cosh⁡⁡ζscosh⁡⁡ζs-cos⁡⁡θ)=cm33π(-2sinh⁡⁡ζs∂∂ζs1cosh⁡⁡ζs-cos⁡⁡θ-cosh⁡⁡ζscosh⁡⁡ζs-cos⁡⁡θ)=-2cm33π∑l=0∞(2sinh⁡⁡ζs∂∂ζs+cosh⁡⁡ζs)e-(l+1/2)|ζs|Pl(cos⁡⁡θ)=-2cm33π∑l=0∞(cosh⁡⁡ζs-(2l+1)sinh⁡⁡|ζs|)e-(l+1/2)|ζs|Pl(cos⁡⁡θ)
for κ=3. The constant coefficients aℓ(ie)1,aℓ(ie)2 for ℓ≥1 and aℓ(ie)3 for ℓ≥0 are to be evaluated from the two-set boundary conditions (4.6) and (4.7), applied on ζs≡-ζ1,ζ2. 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 six-term 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(4.8)Ll(ζs)=e((2l+1)/2)ζsforl≥0,ζs≡-ζ1,ζ2,
then in terms of (4.8) on surfaces S1 at ζs≡-ζ1 and S2 at ζs≡ζ2, we introduce(4.9)Dl,Sjκ(ζs)=22cmκsinh⁡⁡ζs3πLl(ζs)forκ=1,2,Dl,Sj3(ζs)=22cm33πLl(ζs)(lsinh⁡⁡|ζs|-e-|ζs|2)(κ=3)
as well as the easy-to-handle functions, independent of κ=1,2,3,(4.10)Bl,Sj(i)(ζs)=-[lLl-1(ζs)+(l+1)Ll+1(ζs)],Bl,Sj(e)(ζs)=lLl-1(ζs)+l+1Ll+1(ζs)forl≥0,
followed by the simple notations(4.11)Al,Sj(i)κ(ζs)=lLl(ζs)forκ=1,2,Al,Sj(i)3(ζs)=(l+1)Ll(ζs)(κ=3),Cl,Sj(i)κ(ζs)=(l+1)Ll(ζs)forκ=1,2,Cl,Sj(i)3(ζs)=lLl(ζs)(κ=3),
corresponding to interior constants and(4.12)Al,Sj(e)κ(ζs)=-lLl(ζs)forκ=1,2,Al,Sj(e)3(ζs)=-l+1Ll(ζs)(κ=3),Cl,Sj(e)κ(ζs)=-l+1Ll(ζs)forκ=1,2,Cl,Sj(e)3(ζs)=-lLl(ζs)(κ=3),
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 P01(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(4.13)Ãlκxlκ=blκforl=0,1,2,…,κ=1,2,3,
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(4.14)Ãlκ=[B0,S1(i)B0,S1(e)C1,S1(i)κC1,S1(e)κ00000⋯0000000⋯B0,S2(i)B0,S2(e)C1,S2(i)κC1,S2(e)κ00000⋯0000000⋯A0,S1(i)κA0,S1(e)κB1,S1(i)B1,S1(e)C2,S1(i)κC2,S1(e)κ000⋯0000000⋯A0,S2(i)κA0,S2(e)κB1,S2(i)B1,S2(e)C2,S2(i)κC2,S2(e)κ000⋯0000000⋯00A1,S1(i)κA1,S1(e)κB2,S1(i)B2,S1(e)C3,S1(i)κC3,S1(e)κ0⋯0000000⋯00A1,S2(i)κA1,S2(e)κB2,S2(i)B2,S2(e)C3,S2(i)κC3,S2(e)κ0⋯0000000⋯⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮000000000⋯Al-1,S1(i)κAl-1,S1(e)κBl,S1(i)Bl,S1(e)Cl+1,S1(i)κCl+1,S1(e)κ0⋯000000000⋯Al-1,S2(i)κAl-1,S2(e)κBl,S2(i)Bl,S2(e)Cl+1,S2(i)κCl+1,S2(e)κ0⋯⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮]and is of infinite dimensions. The matrix of constant terms bℓκ is given by(4.15)blκ=[D0,S1κD0,S2κD1,S1κD1,S2κ⋯Dl,S1κDl,S2κ⋯]⊺,
and the matrix of unknown constant coefficients xℓκ renders(4.16)xlκ=[a0(i)κa0(e)κa1(i)κa1(e)κ⋯al-1(i)κal-1(e)κal(i)κal(e)κal+1(i)κal+1(e)κ⋯]⊺,
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 (4.15) one have to set D0,S1κ=D0,S2κ≡0, while in (4.16) the constant coefficients a0(i)κ and a0(e)κ do not affect (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ℓ(ie)1,aℓ(ie)2 for ℓ≥1 and aℓ(ie)3 for ℓ≥0 can be obtained up to any order of accuracy through cut-off techniques based on the solution of square systems. Once the solution of the system is obtained, the scattered magnetic field (4.2), defined in the domain (3.3), is provided by(4.17)H3s=∇{∑l=1∞cosh⁡⁡ζ-cos⁡⁡θ×[∑l=1∞(al(i)1e((2l+1)/2)ζ+al(e)1e-((2l+1)/2)ζ)Pl1(cos⁡⁡θ)cos⁡⁡φ+∑l=1∞(al(i)2e((2l+1)/2)ζ+al(e)2e-((2l+1)/2)ζ)Pl1(cos⁡⁡θ)sin⁡⁡φ+∑l=0∞(al(i)3e((2l+1)/2)ζ+al(e)3e-((2l+1)/2)ζ)Pl(cos⁡⁡θ)]}forr∈V(R3)-{r0},
reminding that for n=3 the scattered electric field vanishes, E3s=0.The efficiency of the cut-off method is demonstrated in Section5 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 whenn=0, that is, the static term H0s, though not an easy task as immediate consequence of the complexity of the incident field H0i, given by (2.10). This field admits the double action of the gradient operator (at position r≠r0) on the quantity 1/R for R=|r-r0|. 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 S1 and S2 for n=0, whereas for n̂ defined by (3.6) with (3.7), it is stated by(4.18)H0s=∇Φ0swithΔΦ0s=0,whereζ̂⋅(H0i+H0s)=0on⁡ζ=ζs≡-ζ1,ζ2.
Similarly with the previous analysis, the interior and exterior harmonic potential Φ0s, with the help of relations (3.13)–(3.15), enter into (4.18) to yield(4.19)H0s=∑l=0∞∑m=0l∑q=e,o∇[cosh⁡⁡ζ-cos⁡⁡θ(bl(i)mqe(l+1/2)ζ+bl(e)mqe-(l+1/2)ζ)Plm(cos⁡⁡θ)fmq(φ)]
for every r∈V(ℝ3)-{r0}, where, now, bℓ(i)mq and bℓ(e)mq for ℓ=0,1,2,…,m=0,1,2,…,ℓ, and q=e,o denote the constant coefficients to evaluate by the appropriate boundary condition (4.18). Initially, we calculate the two summands of the condition separately. Then, in view of expression (3.10), we come up with (4.20)ζ̂⋅H0s|ζ=ζs={cosh⁡⁡ζs-cos⁡⁡θ2c×∑l=0∞∑m=0l∑q=e,o[bl(i)mq(sinh⁡⁡ζs+(2l+1)(cosh⁡⁡ζs-cos⁡⁡θ))e(l+1/2)ζs+bl(e)mq(sinh⁡⁡ζs-(2l+1)(cosh⁡⁡ζs-cos⁡⁡θ))e-(l+1/2)ζs]×Plm(cos⁡⁡θ)fmq(φ)∑m=0l}
for ζs≡-ζ1,ζ2. 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 H0s. 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(4.21)ζ̂⋅H0i|ζ=ζs=cosh⁡⁡ζs-cos⁡⁡θc[∂∂ζ(∇r1|r-r0|)]ζ=ζs⋅m4π,
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(4.22)ζ̂⋅H0i|ζ=ζs=-cosh⁡⁡ζs-cos⁡⁡θc[∂∂ζ(∇r01|r-r0|)]ζ=ζs⋅m4π,
which, upon introduction of eigenexpansion (3.20) and with the help of the derivative (3.19), becomes(4.23)ζ̂⋅H0i|ζ=ζs=-cosh⁡⁡ζs-cos⁡⁡θ2c×∑l=0∞∑m=0l∑q=e,o{(m4π⋅∇r0ρlmq(ζs;r0))×[(sinh⁡⁡ζs+sgn⁡(ζ0-ζs)(2l+1)(cosh⁡⁡ζs-cos⁡⁡θ))esgn⁡(ζ0-ζs)(l+1/2)ζs]×Plm(cos⁡⁡θ)fmq(φ)(m4π⋅∇r0ρlmq(ζs;r0))},
where ρℓmq(ζs;r0) and sign function sgn(ζ0-ζs) are provided by relation (3.17) and (3.18), respectively, on surfaces ζ=ζs. The gradient ∇r0ρℓmq(ζs;r0) is a known quantity, which can be determined with the aid of (3.10), (3.18)-(3.19) taken at r=r0, 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 H3s 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 S1 at ζs≡-ζ1 and S2 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(4.24)Dl,Sjmq(ζs;r0)=Fl(ζs;ζ0)Mlmq(ζs;r0)+(l-m)Gl-1(ζs;ζ0)Ml-1mq(ζs;r0)+(l+m+1)Gl+1(ζs;ζ0)Ml+1mq(ζs;r0),
in terms of functions(4.25)Fl(ζs;ζ0)=-[sinh⁡⁡ζs+sgn⁡(ζ0-ζs)(2l+1)cosh⁡⁡ζs]esgn⁡(ζ0-ζs)((2l+1)/2)ζs,Gl(ζs;ζ0)=sgn⁡(ζ0-ζs)esgn⁡(ζ0-ζs)((2l+1)/2)ζs,Mlmq(ζs;r0)=m4π⋅∇r0ρlmq(ζs;r0),where∇r0=cosh⁡⁡ζ0-cos⁡⁡θ0c[ζ̂0∂∂ζ0+θ̂0∂∂θ0+φ̂0sin⁡⁡θ0∂∂φ0]
and ℓ=0,1,2,…,m=0,1,2,…,ℓ, and q=e,o. In addition, we recall definition (4.8) from the previous problem for H3s(4.26)Ll(ζs)=e((2l+1)/2)ζsforl≥0,ζs≡-ζ1,ζ2,
so as to obtain the notations (4.10), which are(4.27)Bl,Sj(i)(ζs)=-[lLl-1(ζs)+(l+1)Ll+1(ζs)],Bl,Sj(e)(ζs)=lLl-1(ζs)+l+1Ll+1(ζs)forl≥0,
as well as the easy-to-handle functions(4.28)Al,Sj(i)m(ζs)=(l-m+1)Ll(ζs),Al,Sj(e)m(ζs)=-l-m+1Ll(ζs)forl=0,1,2,…,m=0,1,2,…,l,Cl,Sj(i)m(ζs)=(l+m)Ll(ζs),Cl,Sj(e)m(ζs)=-l+mLl(ζs)forl=0,1,2,…,m=0,1,2,…,l,
corresponding to interior and exterior constants, respectively. Consequently, the similar matrix form of the systems of the linear algebraic equations discussed earlier is(4.29)Ãlmxlmq=blmqforl=0,1,2,…,m=0,1,2,…,l,q=e,o,
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(4.30)Ãlm=[B0,S1(i)B0,S1(e)C1,S1(i)mC1,S1(e)m00000⋯0000000⋯B0,S2(i)B0,S2(e)C1,S2(i)mC1,S2(e)m00000⋯0000000⋯A0,S1(i)mA0,S1(e)mB1,S1(i)B1,S1(e)C2,S1(i)mC2,S1(e)m000⋯0000000⋯A0,S2(i)mA0,S2(e)mB1,S2(i)B1,S2(e)C2,S2(i)mC2,S2(e)m000⋯0000000⋯00A1,S1(i)mA1,S1(e)mB2,S1(i)B2,S1(e)C3,S1(i)mC3,S1(e)m0⋯0000000⋯00A1,S2(i)mA1,S2(e)mB2,S2(i)B2,S2(e)C3,S2(i)mC3,S2(e)m0⋯0000000⋯⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮000000000⋯Al-1,S1(i)mAl-1,S1(e)mBl,S1(i)Bl,S1(e)Cl+1,S1(i)mCl+1,S1(e)m0⋯000000000⋯Al-1,S2(i)mAl-1,S2(e)mBl,S2(i)Bl,S2(e)Cl+1,S2(i)mCl+1,S2(e)m0⋯⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮]of infinite dimensions. For ℓ≥0,m=0,1,2,…,ℓ and q=e,o the matrix of the constant terms bℓmq is provided by(4.31)blmq=[D0,S1mqD0,S2mqD1,S1mqD1,S2mq⋯Dl,S1mqDl,S2mq⋯]⊺
and the matrix of the unknown constant coefficients xℓmq is assumed to be(4.32)xlmq=[b0(i)mqb0(e)mqb1(i)mqb1(e)mq⋯bl-1(i)mqbl-1(e)mqbl(i)mqbl(e)mqbl+1(i)mqbl+1(e)mq⋯]⊺,
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ℓ(i)mq 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(4.33)H0s=∇[cosh⁡ζ-cos⁡θ∑l=0∞∑m=0l∑q=e,o(bl(i)mqe((2l+1)/2)ζ+bl(e)mqe-((2l+1)/2)ζ)Plm(cos⁡⁡θ)fmq(φ)]
for every r∈V(ℝ3)-{r0}, while for n=0 there exists no scattered electric field, E0s=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 atn=2, where a very difficult and cumbersome manipulation of the boundary value problem (2.20) with (2.23) results into the dynamic scattered fields H2s and E2s.There exist two reasons that are responsible for this difficulty. The first one is the coupling of the particular model with the zero-order fieldH0s (static term), while the second one refers to the extra electric field E2s, which enters our problem with the corresponding additional boundary conditions. However, H2s as well as E2s 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 imaginary-valued (quadrature) part of the magnetic field Hs and to the entire imaginary term of the corresponding electric field Es. Indeed, only the real-valued (in-phase) part of Hs is essentially made of the static contribution H0s. 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.34)H2s=Φ2s+12(rΦ0s)withΔΦ2s=0,ΔΦ0s=0,whereζ̂⋅(H2i+H2s)=0on⁡ζ=ζs≡-ζ1,ζ2,(4.35)E2s=1σ∇×H2s=1σ(∇×Φ2s+12∇Φ0s×r),whereζ̂×(E2i+E2s)=0on⁡ζ=ζs≡-ζ1,ζ2,
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 E2s is obvious, this is not the case for the scattered magnetic field H2s, where we have(4.36)∇⋅H2s=0⟹2∇⋅Φ2s+3Φ0s+(r⋅∇Φ0s)=0,
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Φ0s), it is ensured from identity (A.11) as well as from the harmonic character of position vector r and potential Φ0s that(4.37)Δ(rΦ0s)=2∇Φ0s≡2H0s,
since ∇⊗r=Ĩ, where coupling with the n=0 solution is exhibited for the nonharmonic part of the field H2s. 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 Φ0s as(4.38)Φ0s=cosh⁡⁡ζ-cos⁡⁡θ∑l=0∞∑m=0l∑q=e,o(bl(i)mqe(l+1/2)ζ+bl(e)mqe-(l+1/2)ζ)Plm(cos⁡⁡θ)fmq(φ)
for r∈V(ℝ3)-{r0}. Additionally, the second set of functions of the dynamic field H2s is built up from the harmonic character of Φ2s for internal and external domains, that is,(4.39)Φ2s=cosh⁡⁡ζ-cos⁡⁡θ∑l=0∞∑m=0l∑q=e,o(hl(i)mqe(l+1/2)ζ+hl(e)mqe-(l+1/2)ζ)Plm(cos⁡⁡θ)fmq(φ),
where(4.40)hl(ie)mq=cl(ie)mqx̂1+dl(ie)mqx̂2+el(ie)mqx̂3forl=0,1,2,…,m=0,1,2,…,l,q=e,o
denote the vector character of the unknown constant coefficients hℓ(ie)mq, 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(4.41)r=∑i=13xix̂i=-ccosh⁡⁡ζ-cos⁡⁡θ(sinh⁡⁡ζcos⁡⁡θζ̂+cosh⁡⁡ζcos⁡⁡θθ̂)
for every ζ∈ℝ and θ∈[0,π], the full solution (4.34) for the scattered magnetic field is expressed as(4.42)H2s=cosh⁡⁡ζ-cos⁡⁡θ×∑l=0∞∑m=0l∑q=e,o[(hl(i)mq+bl(i)mq2r)e(l+1/2)ζ+(hl(e)mq+bl(e)mq2r)e-(l+1/2)ζ]×Plm(cos⁡⁡θ)fmq(φ),
while the scattered electric field (4.35) assumes the form(4.43)E2s=1σ∑l=0∞∑m=0l∑q=e,o[∇(cosh⁡⁡ζ-cos⁡⁡θe(l+1/2)ζPlm(cos⁡⁡θ)fmq(φ))×(hl(i)mq+bl(i)mq2r)+∇(cosh⁡⁡ζ-cos⁡⁡θe-(l+1/2)ζPlm(cos⁡⁡θ)fmq(φ))×(hl(e)mq+bl(e)mq2r)],
by straightforward application of identity (A.4), since ∇×r=0. Obviously, both electromagnetic fields for n=2 are defined in the domain V(ℝ3)-{r0} (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ℓ(ie)mq, which for the prescribed ζ∈(-ζ1,ζ2), θ∈[0,π] and φ∈[0,2π) is(4.44)hl(ie)mq=1cosh⁡⁡ζ-cos⁡⁡θ×{[-sinh⁡ζsin⁡θ(cl(ie)mqcos⁡φ+dl(ie)mqsin⁡φ)-el(ie)mq(cosh⁡ζcos⁡θ-1)]ζ̂+[(cosh⁡ζcos⁡θ-1)(cl(ie)mqcos⁡φ+dl(ie)mqsin⁡φ)-el(ie)mqsinh⁡ζsin⁡θ]θ̂+[(cosh⁡ζ-cos⁡θ)(-cl(ie)mqsin⁡φ+dl(ie)mqcos⁡φ)]φ̂},
where ℓ=0,1,2,…,m=0,1,2,…,ℓ, and s=e,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ℓ(ie)mq, dℓ(ie)mq, and eℓ(ie)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Φ0s) 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 forn=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 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 S1:ζ=ζs≡-ζ1 and S2:ζ=ζ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:(4.45)∑l=0∞∑m=0l∑q=e,o[fl,c(i)mq,κ(ζs,θ,φ)cl(i)mq+fl,d(i)mq,κ(ζs,θ,φ)dl(i)mq+fl,e(i)mq,κ(ζs,θ,φ)el(i)mq+fl,c(e)mq,κ(ζs,θ,φ)cl(e)mq+fl,d(e)mq,κ(ζs,θ,φ)dl(e)mq+fl,e(e)mq,κ(ζs,θ,φ)el(e)mq]=∑l=0∞∑m=0l∑q=e,oglmq,κ(ζs,θ,φ;r0)forκ=1,2,3,ζs≡-ζ1,ζ2,
where θ∈[0,π],φ∈[0,2π), and r0=(ζ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:(4.46)fl,c(i)mq,1(ζs,θ,φ)=2sinh⁡ζse((2l+1)/2)ζssin⁡θcos⁡φYlmq(θ,φ),fl,c(e)mq,1(ζs,θ,φ)=e-(2l+1)ζsfl,c(i)mq,1(ζs,θ,φ),fl,d(i)mq,1(ζs,θ,φ)=2sinh⁡ζse((2l+1)/2)ζssin⁡θsin⁡φYlmq(θ,φ),fl,d(e)mq,1(ζs,θ,φ)=e-(2l+1)ζsfl,d(i)mq,1(ζs,θ,φ),fl,e(i)mq,1(ζs,θ,φ)=2(cosh⁡ζscos⁡θ-1)e((2l+1)/2)ζsYlmq(θ,φ),fl,e(e)mq,1(ζs,θ,φ)=e-(2l+1)ζsfl,e(i)mq,1(ζs,θ,φ),glmq,1(ζs,θ,φ;r0)=(cosh⁡ζs-cos⁡θ)Sl(ζs,ζ0)ζ̂(ζs,θ,φ)⋅[Mlmq(ζs;r0)r0-ρlmq(ζs;r0)m4π]Ylmq(θ,φ)+csinh⁡ζscos⁡θ[Sl(ζs;ζ0)Mlmq(ζs;r0)-Blmq(ζs)]Ylmq(θ,φ),
while for κ=2 it is readily obtained(4.47)fl,c(i)mq,2(ζs,θ,φ)=[2sinh⁡ζssin⁡θPlmq(ζs,θ,φ)e((2l+1)/2)ζs+hl,+(ζs,θ)fl,e(i)mq,1(ζs,θ,φ)]cos⁡φ,fl,c(e)mq,2(ζs,θ,φ)=e-(2l+1)ζs×[2sinh⁡ζssin⁡θPlmq(ζs,θ,φ)e((2l+1)/2)ζs+hl,-(ζs,θ)fl,e(i)mq,1(ζs,θ,φ)]cos⁡φ,fl,d(i)mq,2(ζs,θ,φ)=[2sinh⁡ζssin⁡θPlmq(ζs,θ,φ)e((2l+1)/2)ζs+hl,+(ζs,θ)fl,e(i)mq,1(ζs,θ,φ)]sin⁡φ,fl,d(e)mq,2(ζs,θ,φ)=e-(2l+1)ζs×[2sinh⁡ζssin⁡θPlmq(ζs,θ,φ)e((2l+1)/2)ζs+hl,-(ζs,θ)fl,e(i)mq,1(ζs,θ,φ)]sin⁡φ,fl,e(i)mq,2(ζs,θ,φ)=2e((2l+1)/2)ζs×[(cosh⁡ζscos⁡θ-1)Plmq(ζs,θ,φ)-sinh⁡ζssin⁡θhl,+(ζs,θ)Ylmq(θ,φ)],fl,e(e)mq,2(ζs,θ,φ)=2e-((2l+1)/2)ζs×[(cosh⁡ζscos⁡θ-1)Plmq(ζs,θ,φ)-sinh⁡ζssin⁡θhl,-(ζs,θ)Ylmq(θ,φ)],glmq,2(ζs,θ,φ;r0)=c(bl(i)mqe((2l+1)/2)ζshl,+(ζs,θ)+bl(e)mqe-((2l+1)/2)ζshl,-(ζs,θ))×cosh⁡ζssin⁡θYlmq(θ,φ)-4c(cosh⁡ζs-cos⁡θ)Sl(ζs,ζ0)φ̂(φ)⋅Mlmq(ζs;r0)Ylmq(θ,φ)-cBlmq(ζs)sinh⁡ζscos⁡θPlmq(ζs,θ,φ).
Finally, the κ=3 case is structured by the functions(4.48)fl,c(i)mq,3(ζs,θ,φ)=2sinh⁡ζse((2l+1)/2)ζssin⁡θcos⁡φFlmq(ζs,θ,φ)-cosh⁡ζs-cos⁡θsinh⁡ζshl,+(ζs,θ)fl,d(i)mq,1(ζs,θ,φ),fl,c(e)mq,3(ζs,θ,φ)=e-(2l+1)ζs[-cosh⁡ζs-cos⁡θsinh⁡ζshl,-(ζs,θ)fl,d(i)mq,1(ζs,θ,φ)]2sinh⁡ζse((2l+1)/2)ζssin⁡θcos⁡φFlmq(ζs,θ,φ)-cosh⁡ζs-cos⁡θsinh⁡ζshl,-(ζs,θ)fl,d(i)mq,1(ζs,θ,φ)],fl,d(i)mq,3(ζs,θ,φ)=2sinh⁡ζse((2l+1)/2)ζssin⁡θsin⁡φFlmq(ζs,θ,φ)+cosh⁡ζs-cos⁡θsinh⁡ζshl,+(ζs,θ)fl,c(i)mq,1(ζs,θ,φ),fl,d(e)mq,3(ζs,θ,φ)=e-(2l+1)ζs[-cosh⁡ζs-cos⁡θsinh⁡ζshl,-(ζs,θ)fl,d(i)mq,1(ζs,θ,φ)]2sinh⁡ζse((2l+1)/2)ζssin⁡θsin⁡φFlmq(ζs,θ,φ)+cosh⁡ζs-cos⁡θsinh⁡ζshl,-(ζs,θ)fl,c(i)mq,1(ζs,θ,φ)],fl,e(i)mq,3(ζs,θ,φ)=2(cosh⁡ζscos⁡θ-1)e((2l+1)/2)ζsFlmq(ζs,θ,φ),fl,e(e)mq,3(ζs,θ,φ)=e-(2l+1)ζsfl,e(i)mq,3(ζs,θ,φ),glmq,3(ζs,θ,φ;r0)=4csin⁡θ(cosh⁡ζs-cos⁡θ)Sl(ζs,ζ0)θ̂(ζs,θ,φ)⋅Mlmq(ζs;r0)Ylmq(θ,φ)-cBlmq(ζs)sinh⁡ζscos⁡θFlmq(ζs,θ,φ).
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;r0) 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ℓ(e)mq for ℓ≥0,m=0,1,2,…,ℓ, and q=e,o are taken from the solution of the zero-order problem (4.29)–(4.32). On the other hand, the rest of helpful functions, in the aforementioned formulae, are(4.49)Ylmq(θ,φ)=Plm(cos⁡θ)fmq(φ),Blmq(ζs)=bl(i)mqe((2l+1)/2)ζs+bl(e)mqe-((2l+1)/2)ζs,Mlmq(ζs;r0)=m4π⋅∇r0ρlmq(ζs;r0),Mlmq(ζs;r0)=m4π×∇r0ρlmq(ζs;r0),
where(4.50)∇r0=cosh⁡ζ0-cos⁡θ0c[ζ̂0∂∂ζ0+θ̂0∂∂θ0+φ̂0sin⁡θ0∂∂φ0],hl,±(ζs,θ)=sinh⁡ζs±(2l+1)(cosh⁡ζs-cos⁡θ),Sl(ζs,ζ0)=esgn⁡(ζ0-ζs)((2l+1)/2)ζs,
while(4.51)Plmq(ζs,θ,φ)=sin⁡θ[Plm(cos⁡θ)-2(cosh⁡ζs-cos⁡θ)∂Plm(cos⁡θ)∂(cos⁡θ)]flmq(φ),Flmq(ζs,θ,φ)=2Plm(cos⁡θ)(cosh⁡ζs-cos⁡θ)∂flmq(φ)∂φ.
The form of conditions (4.45) is certainly not appropriate to apply orthogonality rules neither in the θ-dependence via (A.19) nor in the φ-dependence. 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 ∑ℓ=0∞∑m=0ℓ∑q=e,o(const.)Pℓm(cosθ)fmq(φ) 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 cℓ(ie)mq,dℓ(ie)mq, and eℓ(ie)mq, therefore the vector ones hℓ(ie)mq=cℓ(ie)mqx̂1+dℓ(ie)mqx̂2+eℓ(ie)mqx̂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.


## 5. 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 acceptsR-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 zero-order 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 Hs≡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 x3-axis at distance d>α1+α2, while embedded in a homogeneous infinite space of magnetic permeability μ=4π·10-7F/m and of electric conductivity σ=2·10-3S/m. A magnetic dipole as (2.2), that is, m=m3x̂3 of strength m3=4π·103A·m2, is located nearby at (x10,x20,x30) and illuminates the two spheres at the low frequency of 500Hz, while the magnetic field is measured along a line (a borehole) at (x1p,x2p,x3p). 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 Figure2, where in accordance with the position of the dipole, the results are symmetric with respect to the line x3=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 x3>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 real and the imaginary parts of the approximated low-frequency scattered magnetic fieldH for the orders n=0,2,3. The two spheres have equal radii α1=α2=25m and their centers are at d=75m. The magnetic dipole is located at (x10,x20,x30)=(200m,0m,0m) and the measurement line probe is placed at (x1p,x2p,x3p)=(0m,200m,[-400,400]m), obtaining the field on the x3-axis.
(a)(b)(c)(d)(e)(f)The real and the imaginary parts of the approximated low-frequency scattered magnetic fieldH for the orders n=0,2,3. The two spheres have different radii α1=25m and α2=50m, while their centers are at the somewhat far distance d=100m. The magnetic dipole is located at (x10,x20,x30)=(200m,0m,0m) and the measurement line probe is placed at (x1p,x2p,x3p)=(0m,200m,[-400,400]m), obtaining the field on the x3-axis.
(a)(b)(c)(d)(e)(f)The real and the imaginary parts of the approximated low-frequency scattered magnetic fieldH for the orders n=0,2,3. The two spheres have different radii α1=25m and α2=50m, while their centers are very close (d=76m). The magnetic dipole is located at (x10,x20,x30)=(200m,0m,0m) and the measurement line probe is placed at (x1p,x2p,x3p)=(0m,200m,[-400,400]m), obtaining the field on the x3-axis.
(a)(b)(c)(d)(e)(f)The series solution obtained for the evaluation ofH 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 dependson 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.


## 6. 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.



---

*Source: 628261-2012-04-02.xml*