A Numerical Alternative for 3D Addition Theorems Based on the Bilinear Form of the Dyadic Green's Function and the Equivalence Principle
Giacomo Giannetti
Ludger Klinkenbusch
A numerical method based on the equivalence principle and the dyadic Green's function is presented. It can be used to compute the sphericalmultipole amplitudes with respect to an origin in a subdomain 2 due to sources in a distinct subdomain 1. As an example, consider that subdomain 1 contains a horn antenna that is solved numerically using a commercial fullwave simulator. The radiated field serves as the incident field for subdomain 2 which contains the scatterer, in our example a lossless dielectric sphere. The proposed method is based on the equivalence currents on a Huygens surface enclosing the antenna and uses the freespace dyadic Green's function to compute the electric and magnetic fields on a sphere enclosing the scatterer. From this electromagnetic field on the spherical surface, the sphericalmultipole amplitudes of the incident field with respect to the center of the sphere enclosing the scatterer are obtained numerically and can be further processed. The results obtained with this method are compared to the results solely computed by the numerical fullwave simulator.
 Article
(4273 KB)  Fulltext XML
 BibTeX
 EndNote
Computational Electromagnetics plays a crucial role in electrical engineering, with many applications in microwave techniques, antennas, and propagation, among others (Davidson, 2010; Sumithra and Thiripurasundari, 2017). Several methods have been developed to solve electrically large problems, where the typical size of the structure exceeds tens or even hundreds of wavelengths. In this context, the case of objects separated by a homogeneous background medium is of particular interest and will be addressed in the following.
A typical approach to solving such problems is based on a physical domain decomposition, that is, to first solve the electromagnetic problem in each of the subdomains separately and to subsequently find the solution of the entire problem by letting the subdomain solutions interact through a suitable coupling procedure. The electromagnetic problem in each of the subdomains and the coupling can be solved by a specialized numerical method.
An example for a related method can be found in a recent work by Losenicky et al. (2021), where the Method of Moments (MoM) and the Tmatrix approach are combined. In particular, there the MoM is used to solve the radiator, an electric dipole, and the fields are projected on a sphere enclosing the radiator to obtain a multipole representation. This is then used to translate the multipoles and solve the scattering problem in another subdomain. However, this approach requires that the subdomains are enclosed by two nonintersecting spheres; therefore, this method does not allow having the antenna and the scatterer in close proximity. This aspect is more limiting when the aspect ratio of the antenna or the scatterer is large, e.g., an electric dipole which is large in one dimension and small in the other two. Other instances of related approaches are found in Alian and Oraizi (2018, 2019), where the addition theorem is combined with the equivalence principle algorithm (EPA) (Li et al., 2006; Li and Chew, 2007) to solve scattering problems with multiple PEC objects. Again, spheres enclosing the objects which must not intersect are needed to apply the method. Additionally, in Alian and Oraizi (2018, 2019) far fields are analyzed; hence no radiator is included in the analysis.
In the present work, we extend the 2D multipole approach we introduced in Giannetti and Klinkenbusch (2023) to the 3D case. The method is based on the equivalence principle and the freespace dyadic Green's function. Two subdomains are considered: one contains an antenna, and the other a scatterer. The classical approach is based on the addition theorem for vector sphericalmultipole functions (VSMFs) and requires the Huygens surface enclosing the radiator to be a sphere. The method described here allows one to enclose the radiator with a Huygens surface of arbitrary shape. Hence, it is more flexible and the distances between the different subdomains can be reduced. The field radiated by the scatterer can then be computed using a sphericalmultipole expansion centered at subdomain 2 and a suitable solver. The scatterer can have an arbitrary geometry – though the analysis naturally simplifies for a spherical scatterer, e.g. the model of the human head as used in Losenicky et al. (2021).
The paper is organized as follows. The problem is outlined in Sect. 2 while the proposed method is described in Sect. 3. First numerical results and some conclusions are shown in Sects. 4 and 5, respectively.
Figure 1 shows the problem and the used notations. The scatterer is enclosed by a spherical surface of radius R_{s}. In the phasor domain and at a time factor e^{+jωt}, the total field for ${r}^{\prime}\ge {R}_{\mathrm{s}}$ is split into an incident and a scattered part. The multipole expansions of the corresponding electrical fields are given by (Klinkenbusch, 2008)
respectively. Here, $Z=\sqrt{\mathit{\mu}/\mathit{\epsilon}}$ is the intrinsic wave impedance of the medium (vacuum is considered in the following) and ${\sum}_{n,m}={\sum}_{n=\mathrm{1}}^{N}{\sum}_{m=n}^{n}$, ${a}_{n,m}^{\mathrm{i}\left(\mathrm{s}\right)}$ and ${b}_{n,m}^{\mathrm{i}\left(\mathrm{s}\right)}$ represent the multipole coefficients of the incident and scattered electromagnetic fields, respectively. The vector sphericalmultipole functions (VSMF) ${\mathit{N}}_{n,m}^{\left(q\right)}\left(\mathit{r}\right)$ and ${\mathit{M}}_{n,m}^{\left(q\right)}\left(\mathit{r}\right)$ are defined by (Klinkenbusch, 2008)
Here, ${w}_{n}^{\left(q\right)}\left(kr\right)=\frac{\mathrm{1}}{kr}\frac{\mathrm{d}}{\mathrm{d}r}\left(r{z}_{n}^{\left(q\right)}\left(kr\right)\right)$, $\widehat{r}$ is the unit vector in the radial direction, and $k=\mathit{\omega}\sqrt{\mathit{\epsilon}\mathit{\mu}}$ is the wavenumber of the homogeneous medium. The superscripts used here, (q)=(1) and (q)=(2), indicate that the radial dependence is given by spherical Bessel functions of the first kind (${z}_{n}^{\left(\mathrm{1}\right)}={j}_{n}$) or by spherical Hankel functions of the second kind (${z}_{n}^{\left(\mathrm{2}\right)}={h}_{n}^{\left(\mathrm{2}\right)}$), respectively. Note that spherical Bessel functions of the first kind are regular everywhere and must be used to represent regular fields at the origin (r=0), while at the given time factor only Hankel functions of the second kind comply with the radiation condition for r→∞.
The transverse spherical multipole functions (TSMFs) ${\mathit{m}}_{n,m}(\mathit{\theta},\mathit{\varphi})$ and ${\mathit{n}}_{n,m}(\mathit{\theta},\mathit{\varphi})$ are defined as
where $\widehat{\mathit{\theta}}$ and $\widehat{\mathit{\varphi}}$ are the unit vectors along θ and ϕ, respectively, and where the surface spherical harmonics ${Y}_{n,m}(\mathit{\theta},\mathit{\varphi})$ are defined by
Here, ${P}_{n}^{m}\left(\mathrm{cos}\right(\mathit{\theta}\left)\right)$ denotes an associated Legendre function of the first kind. Note that ${Y}_{n,m}(\mathit{\theta},\mathit{\varphi})=(\mathrm{1}{)}^{m}{Y}_{n,m}^{\ast}(\mathit{\theta},\mathit{\varphi})$ holds, with the asterisk indicating the complex conjugate.
In the absence of a scatterer, the scattered field vanishes and the total field is identical to the incident field, which is in this case also valid for r<R_{s}.
In case a scatterer is present, the multipole coefficients of the incident and scattered fields are related by the scattering matrix, which fully characterizes the scatterer. For the simple case of a homogeneous isotropic dielectric sphere with wavenumber k_{s} and intrinsic wave impedance Z_{s}, the scattering matrix is diagonal, and the relations between the multipole coefficients of the scattered and incident fields are found as
Moreover, the field inside the dielectric sphere can be expanded using the sphericalmultipole expansion
with the multipole amplitudes
First, we apply the equivalence principle with a Huygens surface that completely encloses the antenna (Fig. 1). There, the equivalent currents, that completely represent the antenna outside its subdomain, are given by (Balanis, 2012)
where Γ is the Huygens surface, and $\widehat{\mathit{n}}$ the unit vector directed outwards (Fig. 1). The equivalent currents Eq. (11) are calculated from the fields delivered by a fullwave simulator.
The electromagnetic field outside of the Huygens surface can be expressed with respect to the coordinate system of subdomain 2 according to (Li et al., 2006; Alvarez et al., 2007; Quijano et al., 2011; Balanis, 2012)
where the integral operators 𝒦 and ℒ are defined by
Here, I_{eq} is either M_{eq} or J_{eq}, and r^{′} and r^{′′} are the observation and source points, respectively, described in the coordinate system of the scatterer (Fig. 1). In Eq. (14), the ∇^{′} operates on r^{′}, and the scalar free space Green's function is given by
The freespace dyadic Green's function $\stackrel{\mathrm{\u203e}}{\stackrel{\mathrm{\u203e}}{G}}({\mathit{r}}^{\prime},{\mathit{r}}^{\prime \prime})$ in Eq. (15) is
where $\stackrel{\mathrm{\u203e}}{\stackrel{\mathrm{\u203e}}{I}}$ represents the unit dyadic, $\mathit{R}={\mathit{r}}^{\prime}{\mathit{r}}^{\prime \prime}$, $R=\left\mathit{R}\right$, and $\widehat{\mathit{R}}=\mathit{R}/R$ is the unit vector pointing from the source point to the observation point. Additionally, we have
Evaluating Eq. (12) on the spherical surface enclosing the scatterer, we get ${\mathit{E}}^{\mathrm{i}}\left({\mathit{r}}^{\prime}\right){}_{{r}^{\prime}={R}_{\mathrm{s}}}$, and the multipole coefficients of the incident field ${a}_{n,m}^{\mathrm{i}}$ and ${b}_{n,m}^{\mathrm{i}}$ are calculated from their transversal (i.e., nonradial) components through exploiting the orthogonality of the VSMFs and TSMFs (Klinkenbusch, 2008)
Alternatively, the coefficients ${a}_{n,m}^{\mathrm{i}}$ and ${b}_{n,m}^{\mathrm{i}}$ can be found from the radial components of the electric and magnetic fields (Klinkenbusch, 2008):
For comparison, the incident electric field in Eqs. (19) and (20) or the incident electric and magnetic fields in Eqs. (21) and (22) can be obtained directly by a fullwave simulator.
4.1 Settings
In the following, a single frequency f_{0}=2 GHz is fixed, which corresponds to a wavelength in vacuum of λ≈150 mm. However, the proposed method can be extended to a frequency range, by repeatedly applying the method to a set of frequency points, and subsequently to the timedomain by applying an inverse Fourier transform.
The scatterer is a lossless isotropic dielectric sphere with radius R_{s}=30 mm and a relative dielectric permittivity ε_{r}=2.2. A horn antenna, optimized for working at f_{0}, is considered and analyzed by the commercial CST © timedomain solver. In Fig. 2, the antenna and its technical drawings are shown. The return loss at f_{0} in the free space is 20.7 dB. For a graphical representation, a possible position of the scatterer is also depicted as a green circle in Fig. 2a. Additionally, the Huygens surface for the proposed method is also drawn: it is a parallelepiped and is called enclosing box in Fig. 2.
The parallelepiped has its sides orthogonal to either ${i}_{{x}^{\prime}}$, ${i}_{{y}^{\prime}}$, or ${i}_{{z}^{\prime}}$, its center in the antenna reference system is located at $({x}_{\mathrm{a}}=\mathrm{0},{y}_{\mathrm{a}}=\mathrm{0},{z}_{\mathrm{a}}=\mathrm{31})$ mm and the lengths of the x,y, and zsides of the parallelepiped are 325.1 mm = 2.17λ, 270.5 mm = 1.80λ, and 692.0 mm = 4.62λ, respectively.
The matrix representations of the integral operators Eqs. (12), (13), (19), and (20) are evaluated numerically by applying the point matching method (Chew, 1995). The distance between the matching points is $\mathit{\lambda}/\mathrm{20}\approx \mathrm{7.5}$ mm and hence the total number of points over which the fields are exported is 18 322. For the value of the truncation limit N in Eq. (1), the following rule of thumb is applied (Hansen, 1988)
For comparison with the classical approach based on the addition theorem (Stein, 1961) for VSMFs (Alian and Oraizi, 2018, 2019; Losenicky et al., 2021), the enclosing sphere with radius R_{a}=385 mm is also drawn in Fig. 2. Note that the enclosing sphere is more extensive than the enclosing box, thus limiting the minimum distance between the antenna's and scatterer's subdomains.
Two positions of the scatterer with respect to the reference system of the antenna are considered:

P1: ${O}^{\prime}=(\mathrm{15},\mathrm{30},\mathrm{700})$ mm, r_{0}=700.8 mm;

P2: ${O}^{\prime}=(\mathrm{0},\mathrm{0},\mathrm{375})$ mm, r_{0}=375.0 mm;
where r_{0} is the distance between the origins of the two reference systems (Fig. 1). The electrical distances kr_{0} are 29.4 and 15.7 for P1 and P2, respectively. Note that for P1, both the proposed method and the classical one based on the addition theorem for VSMFs can be applied. However, P2 can be analyzed only with the method proposed here, since in this case the sphere enclosing the antenna and the one enclosing the scatterer intersect.
A coupling between the scatterer and the antenna is expected, but this is not yet modeled by the proposed method. Therefore, to obtain CST results that are comparable with those of the proposed method, the CST reference solution is computed using the equivalence principle as follows:

step 1: the subdomain of the antenna is solved and the equivalent currents on the bounding box are exported;

step 2: the equivalent currents are loaded in CST as nearfield sources and the antenna is replaced by the Huygens surface (Fig. 3).
Note that the field inside the Huygens surface in Fig. 3 does not vanish, which means that, as expected, there is an interaction between the antenna and the scatterer.
Once the proposed method is validated and no reference solution is further needed, it is sufficient to solve in CST only the subdomain of the antenna and to export the equivalent currents on the bounding box (step 1 of the aforementioned bullet list). The equivalent currents are then used as the input for the proposed method, which solves the scattering problem in the subdomain of the scatterer.
4.2 Comparisons
The multipole coefficients delivered by the proposed method and CST are now compared. Similarly to Hansen (2012), the relative error for this method (TM) and for the CST results is defined as
where c is either a or b, and (s) is either (i) for the incident field or (in) for the field inside the scatterer. In Eq. (24), n and m range from 1 to N and from −n to n, respectively, as for the double summations in Eq. (1). The weighting factor F_{n} in Eq. (25) is necessary to account for the fact that the evaluation of the series Eq. (1) is convergent with n, but it is not for the series formed by the multipole amplitudes c_{n,m}. The reason for this can be found in the behaviour of the spherical Bessel function j_{n}(kr) which converge very fast for increasing n. Hence, we observe for the limit for n→∞ of the spherical Bessel function j_{n}(x):
which thus makes Eq. (25) a suitable normalization factor in Eq. (24). On the other hand, F_{n} is closer to one for small values of n and it is exactly one for n=1. This means that the normalization factor F_{n} accordingly reduces the impact of the normalized multipole amplitudes for increasing values of n.
As an example, we first consider the multipole expansion in P1, with and without the scatterer, and second the multipole expansion in P2 without the scatterer.
For P1, the multipole coefficients of the incident electromagnetic field ${a}_{n,m}^{\mathrm{i}}$ and $(Z/j){b}_{n,m}^{\mathrm{i}}$ are shown in Fig. 4 (first row) while those of the field inside the scatterer ${a}_{n,m}^{\mathrm{in}}$ and $({Z}_{\mathrm{s}}/j){b}_{n,m}^{\mathrm{in}}$ are shown in Fig. 4 (second row). The results agree well, except for slight discrepancies between the two methods that occur when the magnitude of the coefficients is less than 10^{0} V m^{−1}, i.e., become less relevant. The maximum value of the relative error Eq. (24) for ${a}_{n,m}^{\mathrm{i}}$ (${b}_{n,m}^{\mathrm{i}}$) is −30.5 dB (−28.1 dB) and it is −25.2 dB (−31.3 dB) for ${a}_{n,m}^{\mathrm{in}}$ (${b}_{n,m}^{\mathrm{in}}$).
For P1, electric fields are also qualitatively compared. The electric field of the impinging wave on a sphere with radius R_{s} is depicted in Fig. 5. We observe in both cases a similarity of the results from CST and the proposed method.
For P2, the multipole coefficients of the incident field are shown in Fig. 4 (third row). For this case too, the results from the two methods agree well, apart for multipole amplitudes with a magnitude less than $\mathrm{2}\times {\mathrm{10}}^{\mathrm{1}}$ V m^{−1}. The maximum value of the error Eq. (24) for ${a}_{n,m}^{\mathrm{i}}$ (${b}_{n,m}^{\mathrm{i}}$) is −47.8 dB (−36.0 dB).
In Fig. 4 (third column), the error decreases fast for increasing multipole order for the term n^{n} in the denominator of the factor F_{n}. Due to this, the maximum value of the denominator in Eq. (25) is obtained for n=1 and $m=\pm \mathrm{1}$ in the examples analyzed. In addition, the relatively large values of the errors may derive from the limited accuracy of the CST reference solution. To support this, the maximum values of the errors are lower for the scatterer in P2, that is, for the scatterer closer to the antenna and hence a smaller solution domain in CST. The errors may decrease when solving the radiator with a dedicated solver.
We have introduced a method based on the dyadic Green's function and the equivalence principle to represent an electromagnetic field in a coordinate system different from the original one. The method works even when the classical approach based on VSMF translation formulas fails. This feature also allows one to reduce the distance between the antenna and the scatterer, particularly for elongated antennas. The proposed method has been compared to the numerical results purely obtained from the fullwave simulator CST, showing good agreement.
As the next step, the authors intend to solve the full scattering problem by including the electromagnetic interaction between different subdomains.
The code and the data that support the findings of this study are available from the corresponding author, Giacomo Giannetti, upon reasonable request.
GG was involved in Conceptualization, Formal Analysis, Methodology, Software, Validation, Visualization, Writing – original draft, and Writing – review & editing; LK was responsible for Supervision, Project administration, and Writing – review & editing.
At least one of the (co)authors is a member of the editorial board of Advances in Radio Science. The peerreview process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This article is part of the special issue “Kleinheubacher Berichte 2023”. It is a result of the Kleinheubacher Tagung 2023, Miltenberg, Germany, 26–28 September 2023.
This paper was edited by Romanus DyczijEdlinger and reviewed by two anonymous referees.
Alian, M. and Oraizi, H.: Electromagnetic multiple PEC object scattering using equivalence principle and addition theorem for spherical wave harmonics, IEEE Trans. Antennas Propag., 66, 6233–6243, 2018. a, b, c
Alian, M. and Oraizi, H.: Application of equivalence principle for EM scattering from irregular array of arbitrarily oriented PEC scatterers using both translation and rotation addition theorems, IEEE Trans. Antennas Propag., 67, 3256–3267, 2019. a, b, c
Alvarez, Y., LasHeras, F., and Pino, M. R.: Reconstruction of equivalent currents distribution over arbitrary threedimensional surfaces based on integral equation algorithms, IEEE Trans. Antennas Propag., 55, 3460–3468, 2007. a
Balanis, C. A.: Advanced engineering electromagnetics, John Wiley & Sons, ISBN 9780470589489, 2012. a, b
Chew, W. C.: Waves and Fields in Inhomogeneous Media, Electromagnetic waves, IEEE Press, ISBN 9780198592242, 1995. a
Davidson, D.: Computational Electromagnetics for RF and Microwave Engineering, Cambridge University Press, ISBN 9780521518918, 2010. a
Giannetti, G. and Klinkenbusch, L.: Comparative study of two multipolebased numerical methods for 2D fieldtranslation schemes, in: 2023 Kleinheubach Conference, 26–28 September 2023, Miltenberg, Germany, IEEE, 1–4, https://ieeexplore.ieee.org/document/10296703 (last access: 2 September 2024), 2023. a
Hansen, J. E.: Spherical nearfield antenna measurements, in: IEE electromagnetic waves series, P. Peregrinus, ISBN 9780863411106, 1988. a
Hansen, T. B.: Numerical investigation of the systemmatrix method for higherorder probe correction in spherical nearfield antenna measurements, Int. J. Antennas Propag., 2012, 493705, https://doi.org/10.1155/2012/493705, 2012. a
Klinkenbusch, L.: Brief review of sphericalmultipole analysis in radio science, URSI Radio Sci. Bull., 2008, 5–16, 2008. a, b, c, d
Li, M.K. and Chew, W. C.: Wavefield interaction with complex structures using equivalence principle algorithm, IEEE Trans. Antennas Propag., 55, 130–138, 2007. a
Li, M.K., Chew, W. C., and Jiang, L. J.: A domain decomposition scheme based on equivalence theorem, Microw. Opt. Technol. Lett., 48, 1853–1857, 2006. a, b
Losenicky, V., Jelinek, L., Capek, M., and Gustafsson, M.: Method of moments and Tmatrix hybrid, IEEE Trans. Antennas Propag., 70, 3560–3574, 2021. a, b, c
Quijano, J. L. A., Scialacqua, L., Zackrisson, J., Foged, L. J., Sabbadini, M., and Vecchi, G.: Suppression of undesired radiated fields based on equivalent currents reconstruction from measured data, IEEE Antennas Wireless Propag. Lett., 10, 314–317, 2011. a
Stein, S.: Addition theorems for spherical wave functions, Quart. Appl. Math., 19, 15–24, 1961. a
Sumithra, P. and Thiripurasundari, D.: Review on computational electromagnetics, Adv. Electromagn., 6, 42–55, 2017. a