A Numerical Alternative for 3D Addition Theorems Based on the Bilinear Form of the Dyadic Green's Function and the Equivalence Principle

Adv. Radio Sci., 22, 9–15, 2024 https://doi.org/10.5194/ars-22-9-2024 © Author(s) 2024. This work is distributed under the Creative Commons Attribution 4.0 License. A Numerical Alternative for 3D Addition Theorems Based on the Bilinear Form of the Dyadic Green’s Function and the Equivalence Principle Giacomo Giannetti1,zand Ludger Klinkenbusch2 1Department of Information Engineering, University of Florence, 50139 Florence, Italy 2Institute of Electrical and Information Engineering, Kiel University, 24143 Kiel, Germany zFor this work, Giacomo Giannetti received the Young Scientist Award from the German URSI member committee at the Kleinheubacher Tagung 2023. Correspondence: Giacomo Giannetti ([email protected]) Received: 31 March 2024 – Revised: 30 July 2024 – Accepted: 12 August 2024 – Published: 6 September 2024 Abstract. A numerical method based on the equivalence principle and the dyadic Green’s function is presented. It can be used to compute the spherical-multipole 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 full-wave 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 full-wave simulator. 1 Introduction 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 T-matrix 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 non-intersecting 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 multiPublished by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V. 10 G. Giannetti and L. Klinkenbusch: A Numerical Alternative for 3D Addition Theorems ple 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 free-space 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 spherical-multipole 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 spherical-multipole 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. 2 Formulation of the Problem Figure 1 shows the problem and the used notations. The scatterer is enclosed by a spherical surface of radius Rs. In the phasor domain and at a time factor e+jωt , the total field for r0≥Rsis split into an incident and a scattered part. The multipole expansions of the corresponding electrical fields are given by (Klinkenbusch, 2008) Ei(r0)=X n,m ai n,mN(1) n,m(r0)+Z jX n,m bi n,mM(1) n,m(r0)(1) Es(r0)=X n,m as n,mN(2) n,m(r0)+Z jX n,m bs n,mM(2) n,m(r0), (2) respectively. Here, Z=√µ/ε is the intrinsic wave impedance of the medium (vacuum is considered in the following) and Pn,m =PN n=1Pn m=−n,ai(s) n,m and bi(s) n,m represent the multipole coefficients of the incident and scattered electromagnetic fields, respectively. The vector sphericalmultipole functions (VSMF) N(q) n,m(r)and M(q) n,m(r)are defined by (Klinkenbusch, 2008) M(q) n,m(r)=z(q) n(kr)mn,m(θ,φ) (3) N(q) n,m(r)=− z(q) n(kr) kr n(n +1)Yn,m(θ,φ)ˆr +w(q) n(kr)nn,m(θ,φ). (4) Figure 1. Definition of the problem: subdomain 1 (antenna), subdomain 2 (scatterer), and corresponding notations. Here, w(q) n(kr) =− 1 kr d drrz(q) n(kr),ˆris the unit vector in the radial direction, and k=ω√εµ 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(1) n=jn) or by spherical Hankel functions of the second kind (z(2) n=h(2) n), 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) mn,m(θ,φ) and nn,m(θ,φ) are defined as mn,m(θ,φ) =− 1 sin(θ) ∂Yn,m(θ,φ) ∂φ ˆ θ+∂Yn,m(θ,φ) ∂θ ˆ φ(5) nn,m(θ,φ) =∂Yn,m(θ,φ) ∂θ ˆ θ+1 sin(θ) ∂Yn,m(θ,φ) ∂φ ˆ φ, (6) where ˆ θand ˆ φare the unit vectors along θand φ, respectively, and where the surface spherical harmonics Yn,m(θ,φ) are defined by Yn,m(θ,φ) =s(n −m)! (n +m)! 2n+1 4πPm n(cos(θ))ejmφ.(7) Here, Pm n(cos(θ)) denotes an associated Legendre function of the first kind. Note that Yn,−m(θ,φ) =(−1)mY∗ n,m(θ,φ) 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 < Rs. 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 ksand intrinsic wave impedance Zs, the scattering matrix is diagonal, and the relations between the multipole coefficients of the scattered and incident fields are found Adv. Radio Sci., 22, 9–15, 2024 https://doi.org/10.5194/ars-22-9-2024 G. Giannetti and L. Klinkenbusch: A Numerical Alternative for 3D Addition Theorems 11 as as n,m ai n,m =−bs n,m bi n,m =w(1) n(kRs)jn(ksRs)−Z Zsw(1) n(ksRs)jn(kRs) w(2) n(kRs)jn(ksRs)−Z Zsw(1) n(ksRs)jn(kRs) .(8) Moreover, the field inside the dielectric sphere can be expanded using the spherical-multipole expansion Ein(r0)=X n,m ain n,mN(1) n,m(r0)+Zs jX n,m bin n,mM(1) n,m(r0), (9) with the multipole amplitudes ain n,m =Zs Z ai n,mjn(kRs)+as n,mh(2) n(kRs) jn(ksRs) bin n,m =bi n,mw(1) n(kRs)+bs n,mw(2) n(kRs) w(1) n(ksRs) .(10) 3 Proposed method 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) Jeq =ˆ n×Ha(0), Meq =−ˆ n×Ea(0), (11) where 0is the Huygens surface, and ˆ nthe unit vector directed outwards (Fig. 1). The equivalent currents Eq. (11) are calculated from the fields delivered by a full-wave 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) E(r0)=−jkZLJeq(r00)−KMeq(r00)(12) H(r0)=−jkY LMeq(r00)+KJeq(r00),(13) where the integral operators Kand Lare defined by KIeq(r00)=Z0 ∇0g(r0,r00)×Ieq(r00)dS00 (14) LIeq(r00)=Z0 Gr0,r00·Ieq r00dS00.(15) Here, Ieq is either Meq or Jeq, and r0and r00 are the observation and source points, respectively, described in the coordinate system of the scatterer (Fig. 1). In Eq. (14), the ∇0 operates on r0, and the scalar free space Green’s function is given by gr0,r00=e−jk|r0−r00| 4π|r0−r00|.(16) The free-space dyadic Green’s function G(r0,r00)in Eq. (15) is Gr0,r00=I+1 k2∇0∇0g(R) = 3 k2R2+3j kR −1ˆ Rˆ R +1−j kR −1 k2R2Ig(R), (17) where Irepresents the unit dyadic, R=r0−r00,R=|R|, and ˆ R=R/R is the unit vector pointing from the source point to the observation point. Additionally, we have ∇0gr0,r00=−jk −1 Rg(R) ˆ R.(18) Evaluating Eq. (12) on the spherical surface enclosing the scatterer, we get Ei(r0)r0=Rs , and the multipole coefficients of the incident field ai n,m and bi n,m are calculated from their transversal (i.e., non-radial) components through exploiting the orthogonality of the VSMFs and TSMFs (Klinkenbusch, 2008) ai n,m =− 1 n(n +1) 1 1 kRs d dr(rjn(kr))r0=Rs · 2π Z 0 π Z 0 Ei(r0)r0=Rs·n∗ n,m(θ0,ϕ0)sin(θ0)dθ0dϕ0(19) Z jbi n,m =1 n(n +1) 1 jn(kRs)· 2π Z 0 π Z 0 Ei(r0)r0=Rs·m∗ n,m(θ0,ϕ0)sin(θ0)dθ0dϕ0.(20) Alternatively, the coefficients ai n,m and bi n,m can be found from the radial components of the electric and magnetic fields (Klinkenbusch, 2008): ai n,m =− 1 n(n +1) kRs jn(kRs)· 2π Z 0 π Z 0 Ei(r0)r0=Rs· ˆr0Y∗ n,m(θ0,ϕ0)sin(θ0)dθ0dϕ0(21) bi n,m =− 1 n(n +1) kRs jn(kRs)· 2π Z 0 π Z 0 Hi(r0)r0=Rs· ˆr0Y∗ n,m(θ0,ϕ0)sin(θ0)dθ0dϕ0.(22) 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 full-wave simulator. https://doi.org/10.5194/ars-22-9-2024 Adv. Radio Sci., 22, 9–15, 2024 12 G. Giannetti and L. Klinkenbusch: A Numerical Alternative for 3D Addition Theorems 4 Numerical results 4.1 Settings In the following, a single frequency f0=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 time-domain by applying an inverse Fourier transform. The scatterer is a lossless isotropic dielectric sphere with radius Rs=30 mm and a relative dielectric permittivity εr= 2.2. A horn antenna, optimized for working at f0, is considered and analyzed by the commercial CST © time-domain solver. In Fig. 2, the antenna and its technical drawings are shown. The return loss at f0in 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 ix0,iy0, or iz0, its center in the antenna reference system is located at (xa=0,ya=0,za=−31)mm and the lengths of the x-,y-, and z-sides 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 λ/20 ≈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) N=dksRs+10e=12.(23) 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 Ra=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: O0=(15,30,700)mm, r0=700.8 mm; –P2: O0=(0,0,375)mm, r0=375.0 mm; where r0is the distance between the origins of the two reference systems (Fig. 1). The electrical distances kr0are 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 Figure 2. Horn antenna: (a) model in CST; (b) yz-cross section (not in scale); (c) xz-cross section (not in scale). The dimensions of the feeding rectangular waveguide are those of the WR430 standard. In (b) and (c), the taper angles are in degrees, while the other dimensions are in millimeters. Figure 3. Equivalence principle in CST with scatterer in P2: maximum magnitude of the electric field on xa=0 mm (Fig. 2b); HS stands for Huygens’ surface while SC for scatterer. 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 Adv. Radio Sci., 22, 9–15, 2024 https://doi.org/10.5194/ars-22-9-2024 G. Giannetti and L. Klinkenbusch: A Numerical Alternative for 3D Addition Theorems 13 Figure 4. Coefficients of the multipole expansions: P1 or P2 indicates the position of the scatterer; a is for the magnitude, b the phase, and c the error Eq. (24); S, if present, indicates that the scatterer is considered. The legend of (P1a) is the same for all the graphs in the first two columns [without (with) scatterer ˜ Z=Z(˜ Z=Zs)]; the legend of (P1c) is the same for all the graphs in the third column. 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 near-field 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 c n,m =20 ·log10  c(s) TM n,m −c(s) CST n,m ·Fn maxn,m c(s) CST n,m ·Fn (24) Fn=1 √n1 nn (25) where cis either aor b, and (s) is either (i)for the incident field or (in)for the field inside the scatterer. In Eq. (24), n and mrange from 1 to Nand from −nto n, respectively, as for the double summations in Eq. (1). The weighting factor Fnin 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 cn,m. The reason for this can be found in the behaviour of the spherical Bessel function jn(kr) which converge very fast for increasing n. Hence, we observe for the limit for n→∞ of the spherical Bessel function jn(x): jn(x) ≈1 2√xnex 2nn for n→∞ (26) which thus makes Eq. (25) a suitable normalization factor in Eq. (24). On the other hand, Fnis closer to one for small https://doi.org/10.5194/ars-22-9-2024 Adv. Radio Sci., 22, 9–15, 2024 14 G. Giannetti and L. Klinkenbusch: A Numerical Alternative for 3D Addition Theorems Figure 5. Projection of the incident electric field on the surface r0= Rs, 0 ≤ϕ0≤π/2, 0 ≤θ0≤πfor the scatterer in P1: (a) proposed method; (b) CST. values of nand it is exactly one for n=1. This means that the normalization factor Fnaccordingly 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 ai n,m and (Z/j)bi n,m are shown in Fig. 4 (first row) while those of the field inside the scatterer ain n,m and (Zs/j)bin n,m 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 100V m−1, i.e., become less relevant. The maximum value of the relative error Eq. (24) for ai n,m (bi n,m) is −30.5 dB (−28.1 dB) and it is −25.2 dB (−31.3 dB) for ain n,m (bin n,m). For P1, electric fields are also qualitatively compared. The electric field of the impinging wave on a sphere with radius Rsis 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 2 ×10−1V m−1. The maximum value of the error Eq. (24) for ai n,m (bi n,m) is −47.8 dB (−36.0 dB). In Fig. 4 (third column), the error decreases fast for increasing multipole order for the term nnin the denominator of the factor Fn. Due to this, the maximum value of the denominator in Eq. (25) is obtained for n=1 and m=±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. 5 Conclusions 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 full-wave 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. Code and data availability. The code and the data that support the findings of this study are available from the corresponding author, Giacomo Giannetti, upon reasonable request. Author contributions. 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. Competing interests. 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. Disclaimer. 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. Special issue statement. 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. Adv. Radio Sci., 22, 9–15, 2024 https://doi.org/10.5194/ars-22-9-2024 G. Giannetti and L. Klinkenbusch: A Numerical Alternative for 3D Addition Theorems 15 Review statement. This paper was edited by Romanus DyczijEdlinger and reviewed by two anonymous referees. References 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. 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. Alvarez, Y., Las-Heras, F., and Pino, M. R.: Reconstruction of equivalent currents distribution over arbitrary three-dimensional surfaces based on integral equation algorithms, IEEE Trans. Antennas Propag., 55, 3460–3468, 2007. Balanis, C. A.: Advanced engineering electromagnetics, John Wiley & Sons, ISBN 9780470589489, 2012. Chew, W. C.: Waves and Fields in Inhomogeneous Media, Electromagnetic waves, IEEE Press, ISBN 9780198592242, 1995. Davidson, D.: Computational Electromagnetics for RF and Microwave Engineering, Cambridge University Press, ISBN 9780521518918, 2010. Giannetti, G. and Klinkenbusch, L.: Comparative study of two multipole-based numerical methods for 2D field-translation 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. Hansen, J. E.: Spherical near-field antenna measurements, in: IEE electromagnetic waves series, P. Peregrinus, ISBN 9780863411106, 1988. Hansen, T. B.: Numerical investigation of the system-matrix method for higher-order probe correction in spherical near-field antenna measurements, Int. J. Antennas Propag., 2012, 493705, https://doi.org/10.1155/2012/493705, 2012. Klinkenbusch, L.: Brief review of spherical-multipole analysis in radio science, URSI Radio Sci. Bull., 2008, 5–16, 2008. Li, M.-K. and Chew, W. C.: Wave-field interaction with complex structures using equivalence principle algorithm, IEEE Trans. Antennas Propag., 55, 130–138, 2007. 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. Losenicky, V., Jelinek, L., Capek, M., and Gustafsson, M.: Method of moments and T-matrix hybrid, IEEE Trans. Antennas Propag., 70, 3560–3574, 2021. 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. Stein, S.: Addition theorems for spherical wave functions, Quart. Appl. Math., 19, 15–24, 1961. Sumithra, P. and Thiripurasundari, D.: Review on computational electromagnetics, Adv. Electromagn., 6, 42–55, 2017. https://doi.org/10.5194/ars-22-9-2024 Adv. Radio Sci., 22, 9–15, 2024