Theory of perturbation of electric potential by a 3D object made of an anisotropic dielectric material

The extended boundary condition method (EBCM) was formulated for the perturbation of a source electric potential by a 3D object composed of a homogeneous anisotropic dielectric medium whose relative permittivity dyadic is positive definite. The formulation required the application of Green's second identity to the exterior region to deduce the electrostatic counterpart of the Ewald--Oseen extinction theorem. The electric potential inside the object was represented using a basis obtained by implementing an affine bijective transformation of space to the Gauss equation for the electric field. The EBCM yields a transition matrix that depends on the geometry and the composition of the 3D object, but not on the source potential.


Introduction
Dating back more than a hundred years, the Ewald-Oseen extinction theorem [1,2] states that when a 3D object is illuminated by an incident timeharmonic electromagnetic field, the electric and magnetic surface current densities induced on the exterior side of its surface produce an electromagnetic field that cancels the incident field throughout the interior region of the object. This theorem is a cornerstone of research on electromagnetic scattering, and particularly of the extended boundary condition method (EBCM) [3] since its inception in 1965 [4]. Since a bilinear expansion of the free-space dyadic Green function is known [4,5], the EBCM can be used to investigate scattering by an object composed of any linear homogeneous medium for which a basis exists to represent the electromagnetic field therein. This requirement can be satisfied by isotropic dielectric-magnetic mediums, isotropic chiral mediums, and orthorhombic dielectric-magnetic mediums with gyrotropic magnetoelectric properties [6]. In addition, bases have been numerically synthesized for certain gyrotropic dielectric-magnetic mediums [7][8][9]. The EBCM literature continues to expand as time marches on [3,[10][11][12] Historically, the situation has been markedly different in electrostatics. Green's second identity yields an expression [13,14] that was used by Farafonov [15] in 2014 to formulate the EBCM for the perturbation (i.e., "scattering") of a source electric potential (the electrostatic counterpart of an incident time-harmonic electromagnetic field) by a 3D object composed of a homogeneous isotropic medium. In the Farafonov formulation, Green's second identity is applied separately to the exterior and interior regions. The potentials in the two regions are certainly different but, as the Green function for the Poisson equation is applicable to both regions, there is no need to set up the electrostatic counterpart of the Ewald-Oseen extinction theorem. Eigenfunctions of the Laplace equation are used in infinite-series representations of the source potential everywhere, the perturbation potential in the exterior region, as well as the potential induced inside the isotropic dielectric object, since the right side of Eq. (1) has a bilinear expansion [16] in terms of those eigenfunctions. All three series are suitably terminated and a transition matrix is obtained to characterize the perturbation of the source potential by the object. Although the Farafonov formulation does not require the object to be axisymmetric, it has been used only for such objects [17][18][19]. Our objective in this paper is to generalize the EBCM for electrostatic problems in which the perturbing 3D object is composed of a homogeneous anisotropic dielectric medium whose permittivity dyadic is positive definite. We use Green's second identity only in the exterior region and deduce the electrostatic counterpart of the Ewald-Oseen extinction theorem therefrom. As the Laplace and Poisson equations apply in the exterior region, our representations of the source and perturbation potentials are the same as in the Farafonov formulation. But, the Laplace equation is inapplicable in any region occupied by a homogeneous anisotropic dielectric medium. Hence, an affine transformation of space in the Gauss equation for the electric field is implemented in order to generate a basis for representing the electric potential induced inside the object [20]. This basis is premised on the eigenfunctions of the Laplace equation in the spherical coordinate system. There is no requirement of axisymmetry in our formulation, as is explicitly demonstrated with numerical results for perturbation by ellipsoids.
The plan of this paper is as follows. Section 2 describes the boundary-value problem, Sec. 3 contains the derivation of integral equations underlying the EBCM for electrostatics, Sec. 4 describes the formulation of the transition matrix that relates the expansion coefficients of the perturbation potential to those of the source potential, and Sec. 5 presents and discusses illustrative numerical results. The paper concludes in Sec. 6 with remarks pertaining to future work.

Boundary-Value Problem
Let the region V in be the interior of the surface S, whereas the region V out be bounded by the surfaces S and S ∞ , as shown in Fig. 1. The coordinate system has its origin lying inside V in , the sphere with the origin as its center and inscribed in V in has radius r in , the sphere circumscribing V in has radius r out , and S ∞ is a sphere of radius r ∞ . The surface S is described by a continuous and once-differentiable function.
The electric potential Φ(r) satisfies the Poisson equation in the region V out , with ε 0 as the permittivity of free space (i.e., vacuum) and the charge density ρ(r) being non-zero only for r ∈ V so ⊂ V out . Every point in V so is at least a distance r so away from the origin, as shown in Fig. 1.
The charge-free region V in is filled with an anisotropic dielectric medium with permittivity dyadic wherein the scalar ε r > 0 and the diagonal dyadic contains the anisotropy parameters α x > 0 and α y > 0. The constitutive principal axes of this medium are thus parallel tox,ŷ, andẑ. Any real symmetric dyadic can be written in the form of Eqs. (3) and (4) by virtue of the principal axis theorem [21,22]. Anisotropic materials exist in nature [23], and homogenizable composite materials of this kind can be engineered by properly dispersing dielectric fibers in some host isotropic dielectric material [24].

Integral Equations
Application of Green's second identity to the exterior region V out yields [13,25] where Φ + (r ) is the value of Φ(r ) on the exterior side of S, δ(r − r ) is the Dirac delta, and the unit normal vectorn(r) at r ∈ S ∪ S ∞ points into V out .  Figure 1: Schematic of the boundary-value problem. The region V out is vacuous and the region V in is filled with an anisotropic dielectric material of relative permittivity dyadic ε rel , given by Eqs. (3) and (4).
For finite r and with S ∞ taken to be a spherical surface of radius r ∞ , the potential Φ(r ) at r ∈ S ∞ drops off as 1/r as r ∞ → ∞. Accordingly, the integral on S ∞ in Eq. (5) vanishes. Furthermore, let us define the source potential because it exists everywhere when V in is vacuous (just like V out actually is). Equation (5) then delivers the twin equations and Equation (7a) is the electrostatic equivalent of the Ewald-Oseen extinction theorem [1][2][3], and is the cornerstone of our EBCM formulation. It is convenient to rewrite it as According to Eq. (7b) the electric potential in V out has two components, one of which is the source potential Φ so (r) and the other is the perturbation potential Φ pert (r) due to the region V in being non-vacuous [13] With the usual requirement of the electric field having a finite magnitude everywhere on S, the potential must be continuous across that interface; hence, where Φ − (r ) is the value of Φ(r ) on the interior side of S. Likewise, with the assumption of S being charge-free, the normal component of the electric displacement must be continuous across that surface; hence, Substitution of Eqs. (9a) and (9b) in Eqs. (8a) and (8b) delivers Since Φ − (r) at r ∈ S is nothing but the internal potential Φ int (r) evaluated on the interior side of S, we finally get Knowing Φ so , we first have to find Φ int everywhere in V in after choosing r ∈ V in in Eqs. (11). Thereafter, knowing Φ int , we can find Φ pert everywhere in V out after choosing r ∈ V out in Eqs. (11). Thus, the electrostatic counterpart of the Ewald-Oseen formulation is a crucial ingredient in our EBCM formulation, which follows the style of the Waterman formulation for electromagnetics [4,5].

Source and perturbation potentials
The solutions of the Laplace equation in the spherical coordinate system being known [16], the source potential can be represented as Here, the expansion coefficients A smn are supposed to be known, is a normalization factor with δ mm as the Kronecker delta, and the spherical harmonics involve the associated Legendre function P m n (cos θ) [16]. The definitions of the spherical harmonics mandate that A o0n = 0 ∀n ∈ [0, ∞). Likewise, the perturbation potential can be represented as (15) with unknown expansion coefficients B smn , with B o0n = 0 ∀n ∈ [0, ∞).

Algebraic equations
The Green function G(r, r ) has the bilinear expansion [16] Now, we use the orthogonality relationships of the spherical harmonics on spherical surfaces. First, let the surface r = a in < r in be chosen for applying Eq. (11) top . The source potential on the left side of this equation can be replaced by the series on the right side of Eq. (12). The bilinear expansion of G(r, r ) for r < r has to be used on the right side of Eq. (11) top , since r = a in and r ∈ S. After multiplying both sides of the resulting equation by Y s m n (θ, φ) and integrating over the surface r = a in , we obtain Second, we choose the surface r = a out > r out for applying Eq. (11) bot . The perturbation potential on the left side of this equation can be replaced by the series on the right side of Eq. (15). On the right side of Eq. (11) bot the bilinear expansion of G(r, r ) for r > r must be used, since r = a out and r ∈ S.
Then, after multiplying both sides of the resulting equation by Y s m n (θ, φ) and integrating over the surface r = a out , we get The integral equations (17a) and (17b) incorporate the boundary conditions (9a) and (9b) prevailing on the surface S of the perturbing object V in , but the orthogonalities of the spherical harmonics were not applied on that surface. This is the reason for "extended boundary condition" in the name of EBCM [3,4].

Internal potential
The internal potential Φ int (r) satisfies the equation which emerges from the Gauss equation for the electric field. Equation (18) is not the Laplace equation, although it does simplify to the Laplace equation when α x = α y = 1 (i.e., for an isotropic medium [15,[17][18][19]). After applying a bijective affine transformation of space that maps a sphere into an ellipsoid, the solution of Eq. (18) can be written as the series [20] where Z emn (r) = |A −1 • r| n P m Since the basis formed by the functions Z smn (r) is constructed by a spatial transformation of the eigenfunctions of the Laplace equation in the spherical coordinate system, care must be taken to ensure that the angle lies in the same quadrant as θ, and the angle lies in the same quadrant as φ.

Transition matrix
Substitution of Eq. (19) in Eq. (17a) delivers the algebraic equation ss mm nn C s m n , where the surface integral Likewise, substitution of Eq. (19) in Eq. (17b) leads to ss mm nn C s m n , where the surface integral After the indexes n and n in the foregoing equations are restricted to the range [0, N ], Eq. (22) leads to the matrix equation (in abbreviated notation) whereas Eq. (24) yields the matrix equation Then, where the transition matrix relates the expansion coefficients of the perturbation potential to those of the source potential. Very importantly, this matrix does not depend on the source potential.

Asymptotic expression for perturbation potential
The perturbation potential defined in Eq. (15) can be written as Thus, far away from the object, the perturbation potential can be asymptotically stated as and it is determined only by the four perturbational-potential coefficients B e00 , B e01 , B e11 , and B o11 . We have explicitly retained the two lowest-order terms in Eq. (31), because f (1) pert is not always nonzero. For instance, B e00 = 0 when the source potential varies linearly in a fixed direction and object is axisymmetric as well as composed of an isotropic medium [15,17,19]. Likewise, B e00 = 0 when the object is a sphere and the source is either a point charge or a point dipole, so that f (1) pert = 0 and Φ pert (r) ∝ r −2 as r → ∞ [20]. More generally, the distance after which f (1) pert dominates f (2) pert can depend on the direction, because f (1) pert is independent of θ and φ whereas f (2) pert is not. Indeed, the first term on the right side of Eq. (30c) does not exist for θ = π/2, whereas the second and third terms on the right side of the same equation are absent for θ ∈ {0, π}.

Numerical Results and Discussion
By virtue of Eqs. (23), (25), (28), and (29), the perturbation potential must depend on (i) the geometry of the perturbing object described via S; (ii) the composition of the object, quantitated by ε rel ; and (iii) the spatial profile of the source potential Φ so . We present numerical results in this section to illustrate the effects of S, ε rel , and Φ so on the perturbation potential Φ pert .

Geometry
We chose S to be an ellipsoidal surface defined by the position vector The rotation dyadic is a product of three rotation dyadics, with and Sequentially, R z (α s ) represents a rotation by α s ∈ [0, π] about the z axis, R y (β s ) represents a rotation by β s ∈ [0, π] about the new y axis, and R z (γ s ) represents a rotation by γ s ∈ [0, π] around the new z axis. The shape dyadic is given by Thus, the shape principal axes of the ellipsoidal object are parallel to the unit vectors S •x , S •ŷ , and S •ẑ . The semi-axes of this object are given by aµ parallel to S •x , aν parallel to S •ŷ , and a parallel to S •ẑ , where µ > 0 and ν > 0 are the two aspect ratios of the ellipsoid; thus, r out = max {a, aµ, aν}. Finally, with the geometric mean a ave of the three semi-major axes fixed, the length a = a ave (µν) becomes a function of the aspect ratios µ and ν. Note that a = a ave when µ = ν = 1 (i.e., the ellipsoid reduces to a sphere).

Composition
With the arithmetic mean ε ave of the three eigenvalues of ε rel fixed, the relative permittivity scalar is a function of the anisotropy parameters α x and α y . Note that ε r = ε ave when α x = α y = 1 (i.e., the material becomes isotropic).

Source
Sources of two types were considered for illustrative numerical results: (i) a point-charge source and (ii) a point-dipole source. The coefficient [13] A smn = Q ε 0 for a point charge Q located at and φ o ∈ [0, 2π). For a point dipole p = pp located at the same point, the coefficient [26] A smn = p ε 0 where p > 0; the unit vector p = (x cos φ p +ŷ sin φ p ) sin θ p +ẑ cos θ p contains the angles θ p ∈ [0, π] and φ p ∈ [0, 2π) that define the orientation of the dipole, and ∇ o (...) denotes the gradient with respect to r o .

Convergence
A Mathematica™ program was written to calculate the transition matrix T of an anisotropic ellipsoid for a chosen N . The column vector B containing the perturbation-potential coefficients was then calculated using Eq. (28). Also determined was the spatial profile of Φ pert (r, θ, φ). A convergence test was carried out with respect to N , by calculating the integral at diverse values of r ∈ [1, 10] r out as N was incremented by unity. The iterative process of increasing N was terminated when I(r) converged within a preset tolerance of 1% for a specific value of r. The adequate value of N was higher for lower r, with N = 7 sufficient for r ≥ 1.1r out .

Code validation
Validation of the EBCM code was done in two steps. First, the results for the perturbation of the potential of a point charge by a prolate spheroid (i.e., µ = ν < 1) composed of an isotropic dielectric medium were compared with the corresponding series solution obtained using the eigenfunctions of the Laplace equation in the prolate-spheroidal coordinate system [27]. Figure 2 shows plots of the perturbation potential Φ pert (r, θ, 0) evaluated at r = r/r out ∈ {1.1, 2, 4, 10}. These plots were obtained using the series solution [27] and our EBCM when the source is a point charge Q = 10 −9 /36π C located at r o = 2r out and θ o = 0 (with φ o being irrelevant as sin θ o = 0), whereas the perturbing object is characterized by a ave = 3.82 cm, µ = ν = 2/3, ε ave = 3, and α x = α y = 1. Good agreement can be seen in Fig. 2(a) between the two solutions whenr = 1.1, except in the neighborhood of θ = 0. The difference arises because the spherical harmonics used for representing Φ pert and Φ int do not follow the shape of the prolate spheroid well away from the equator, an issue associated with EBCM even for time-harmonic problems [28,29]. However, the difference diminishes asr increases. Indeed, the difference is vanishingly small forr = 10 in Fig. 2(d). Second, the same comparison was done for a sphere (i.e., µ = ν = 1) made of an anisotropic dielectric medium, the series solution for a sphere being available using the eigenfunctions of the Laplace equation in the spherical coordinate system [20]. The difference in values of Φ pert (r) calculated using the two methods for any r ≥ a was always less than 0.001%.

Numerical results for anisotropic dielectric ellipsoids
Now we present the perturbation potential's variations with respect tor and S when the source is either a point charge or a point dipole. For definiteness, we fixed a ave = 5 cm, µ = 0.8, ν = 1.2, ε ave = 3, α x = 0.5, and α y = 1.5. Also, the point source was located at r o = 2r out , θ o = π/4, and φ o = π/6.

Point-charge source
In order to focus on the effect ofr exclusively, we set S = I (i.e., the shape principal axes coincide with the constitutive principal axes). Figures 3(a-d) present the angular profiles of Φ pert (r r out , θ, φ) forr ∈ {1.1, 2, 4, 10} for a pointcharge source with Q = 10 −10 C, when α s = β s = γ s = 0. The perturbation potential Φ pert (r r out , θ, φ) decreases with increase ofr, as becomes evident by comparing Figs. 3(a), 3(b), 3(c), and 3(d). The angular profiles vary significantly forr ≤ 2, but they change very little asr − 2 increases. We have observed through diverse computations (results not shown) that the foregoing observation is valid even when r o exceeds 2r out . Next, we considered the effect of S = I (i.e., the shape principal axes are rotated with respect to the constitutive principal axes). We repeated the calculations of the previous case when α s = 2π/3, β s = 3π/4, and γ s = 5π/9; the results are depicted in Fig. 4. By comparing Figs. 3 with 4, it can be seen that the fact S = I alters the rate of increase/ decrease of Φ pert (r, θ, φ) in any specific direction. Furthermore, the locations of the extremums of Φ pert (r, θ, φ) in the θφ-plane are affected by the difference S − I whenr < 2; see Figs. 3(a) and 4(a), for instance. The impact of S − I diminishes asr − 2 increases, as can be seen on comparing Figs. 3(c) and 4(c), for instance. Figure 4: Same as Fig. 3, except that α s = 2π/3, β s = 3π/4, and γ s = 5π/9. Figures 5(a-d) present the angular profiles of Φ pert (r r out , θ, φ) forr ∈ {1.1, 2, 4, 10} for a point-dipole source with p = 10 −10 C m, θ p = π/4, and φ p = π/3, when S = I. Just as with the point-charge source in Figs. 3(a-d), Φ pert (r r out , θ, φ) decreases with increase ofr; however, the rate of decrease for the point-dipole source is smaller compared to that for the point-charge source. The effect of the type of source is visually evident on comparing Figs. 3 and 5 for the locations of the extremums of Φ pert (r, θ, φ) in the θφ-plane.

Concluding Remarks
We have formulated the EBCM for the perturbation of a source electric potential by a 3D object composed of a homogeneous anisotropic dielectric medium whose relative permittivity dyadic is positive definite. The electrostatic counterpart of the Ewald-Oseen extinction theorem was derived for that purpose, and the electric potential inside the object was represented using a basis obtained by implementing an affine bijective transformation of space to the Gauss equation for the electric field. This formulation is different from the Farafonov formulation [15,[17][18][19], which is applicable only when the object is composed of a homogeneous isotropic medium.
Even though the object is nonspherical, all potentials were represented in terms of the eigenfunctions of the Laplace equation in the spherical coordinate system. Numerical results have shown that our EBCM formulation fails to yield convergent results when the shape of the perturbing object deviates too much from spherical (i.e., when either |µ−1| is substantial and/or |ν−1| is substantial). The same problem bedevils the analogous formulation of the EBCM even for time-harmonic problems [28,29]. As the Laplace equation is either separable or R-separable in several coordinate systems [30], we plan to use bases in future work that will conform better to the shape of the object.