Inverse Source Solutions With Huygens’ Surface Conforming Distributed Directive Spherical Harmonics Expansions

Time-harmonic inverse source solutions are commonly working with electric and magnetic surface current densities defined and discretized on appropriately chosen Huygens’ surfaces. An efficient meshless alternative are spherical harmonics expansions of low-order distributed along the chosen Huygens’ surface, which still possess pretty good spatial localization properties. Similar to expansions with electric and magnetic surface current densities, distributed expansions with Cartesian spherical harmonics (SHCs) or standard vector spherical harmonics (SHVs) are redundant when they are placed on Huygens’ surfaces. This redundancy is reduced by allowing only spherical harmonics, which can be excited by surface current density distributions in a specific plane, and by forming directive spherical harmonics. This results in a considerable reduction of the necessary number of unknowns for representing the sources and improved conditioning of the discretized integral equations. Inspired by the Huygens’ radiator concept, appropriate directive harmonics on the basis of SHCs and SHVs are derived, implemented, and compared in terms of their performances. The validity of the presented techniques is confirmed via inverse source solutions for synthetic and real measurement data.


I. INTRODUCTION
I NVERSE source solvers have become established tools for antenna field transformations, especially in near-field (NF) antenna measurements, but are, of course, also able to handle far-field (FF) data [1], [2]. In contrast to modal expansion methods [3], inverse source solvers can easily handle irregularly sampled data and arbitrary measurement probes [2], [4]. Moreover, they can consider detailed information about the geometric extent and shape of the antenna under test (AUT), which leads to the particular capability of delivering detailed diagnostic information directly on the geometry model of the AUT [2], [5]. Despite their flexibility and versatility, which is generally accompanied by an increased numerical effort, Manuscript  inverse source solvers can be very efficient when they are accelerated based on the concepts of the fast mulitpole method (FMM) [6], [7] or the multilevel FMM (MLFMM) [8], [9], [10]. Inverse source solvers commonly work with a Huygens' surface resembling the shape of the AUT on which surface current densities are, e.g., represented by Rao-Wilton-Glisson (RWG) basis functions [8]. In MLFMM-accelerated inverse source solvers, which are in the focus of this article, the handling of these basis functions requires substantial computational effort, since the FF radiation patterns of the basis functions need to be precomputed and stored in memory.
During the iterative solution process, these patterns have to be aggregated for all basis functions in a source box on the finest level of the MLFMM octree. These steps can be avoided completely by directly expanding the FF patterns of the source boxes in, e.g., spectral samples or spherical harmonics [1]. Spectral samples are only a good choice if they are combined with an exact global interpolation scheme [10] or if they are directly used to represent the FF pattern of the entire AUT [4]. Scalar spherical harmonics have been used in [1] for the expansion of the Cartesian components (SHCs) of the magnetic vector potential, but transverse electric (TE) and transverse magnetic (TM) vector spherical harmonics (SHVs) as found as Eigen solutions of Maxwell's equations [11] can, of course, be utilized, too. Such expansions of the FF radiation patterns of MLFMM boxes on the finest level will be called distributed spherical harmonics (DSHs) expansions in the following. They are no longer directly connected to a mesh, but a mesh describing the AUT geometry can, of course, still be used to identify those boxes, which are needed to represent the radiation fields of the considered AUT.
For expansions with surface current densities, it is well known from the Huygens' and uniqueness principles [12] that representations with electric and magnetic surface current densities over a Huygens' surface are redundant. Since expansions working only with electric or only with magnetic surface current densities are in general not well conditioned, a popular choice is to impose an additional constraint on the electric and magnetic surface current densities, as, e.g., the Love or zero-field condition [5], [13], [14] or a combinedsource (impedance boundary) condition [15], [16], to reduce the degrees of freedom (DOFs). Similarly, DSH expansions over a Huygens' surface are, of course, also redundant. Scalar spherical harmonics expansions of the Cartesian magnetic 0018-926X © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
vector potential components (SHCs) support all, e.g., three Cartesian electric current density components (not only two as needed) and the SHVs comprising TE-and TM-modes form a complete set of field modes in a similar way as an expansion with electric and magnetic surface current densities. If we only consider those SHs, which can be excited by surface current densities on an assumed Huygens' surface, the number of DOFs needed for a DSH expansion can be reduced considerably. If we need a correct field representation only at "one" side of the Huygens' surface, i.e., in the solution domain, but not in the excluded domain, the resulting set of SHs can be combined with previously removed SH functions to achieve predominant radiation into the solution domain, in a similar way as, e.g., done in [15] for electric and magnetic surface current densities. The goal of this article is, therefore, to improve the efficiency and robustness of DSH expansions, as compared to surface current density expansions, even further by removing the described redundancies of the expansions and by possibly improving the conditioning by directive shaping. To avoid the considerable numerical effort of enforcing a zero-field condition, we transfer the concept of the combined source condition to both SHC and SHV expansions. In both cases, the combined source condition is imposed in a way that basis functions are constructed which exhibit predominant radiation into the solution domain but with reduced radiation toward the excluded domain on the other side of the Huygens' surface. The resulting expansions are called directive SHC (SHC-D) and directive SHV (SHV-D) expansions. Compared to surface current density expansions, the resulting DSHs need in general less DOFs and save, thus, a considerable amount of computational effort and memory, but they have also slightly worse localization capabilities, since there is no specifically defined Huygens' surface. Some first investigations on DSH expansions and their possible improvement have been presented in [17].
In Section II, some defining equations of the inverse source formulation with equivalent surface current densities in spatial and spectral (MLFMM) representation are given, before the standard SHC and SHV expansions are introduced. Based on this, the SHC-D and SHV-D expansions are derived in Section III by imposing the relevant combined source conditions and by identifying the relevant equivalent source variations. In Section IV, the properties and the performances of the new expansions are evaluated and demonstrated for a variety of radiation scenarios, where synthetic and measured NF data is considered. Some conclusions are drawn in Section V.

II. NF TRANSMISSION EQUATION WITH EQUIVALENT CURRENTS AND DSH EXPANSION
The observed probe signal in the considered inverse source scenario as depicted in Fig. 1 is formulated as a spatially weighted radiation integral in the form [1], [2]  where all fields and sources are assumed to follow the suppressed time harmonic dependence e jωt with an angular frequency ω. The equivalent electric and magnetic surface current densities J A and M A represent the radiation behavior of the AUT and reside on a Huygens' surface A, which, in the ideal case, encloses the AUT volume V AUT . The dyadic freespace Green's functions for the electric and magnetic current densities, respectively, G E J and G E M are, e.g., found in the given references.
The observation and source locations are denoted by r and r . The mth measurement reference location is r m , the probe volume is denoted by V w , and w is the spatial probe vector weighting function, representing the receive behavior of the probe.
In a numerical solution of the linear inverse source problem according to the forward operator (1), the spatial surface current densities are commonly discretized according to where β p/q are RWG basis functions [8] and J p and M q are the source unknowns (or DOFs) within the arising linear system of equations. To arrive at an efficient evaluation of the forward operator in (1), we work with the propagating plane-wave spectral formulation [1], [2] which is in our case evaluated in a hierarchical multilevel fashion according to the principles of the MLFMM [10]. The vector k is the wavevector with length k = ω(εμ) 1/2 and unit vectork. The expressionJ +M represents the overall propagating plane-wave spectrum of the AUT computed with respect to a source cluster location r s , comprising the radiation of electric and magnetic surface current densities. The propagating plane-wave expansionw of the spatial probe weighting function w is nothing else than the probe receive pattern given with respect to the observation location r m . The FMM translation operator T L with multipole order L, chosen to achieve a desired accuracy, converts the radiated plane waves from a certain source region into incident plane waves in a distant observation region and is found in the given references.
In the numerical solution of the inverse source problem, we work with an iterative solver and evaluate the forward operator and its adjoint, as needed within the involved matrixvector products, on the fly [2], [18]. As such, the solution domain is subdivided into a regular grid of boxes, hierarchically arranged into an octree according to the principles of the MLFMM. We work here with a bottom-up strategy and fix the box sizes on the finest level, e.g., to a sidelength of half a wavelength. The spectral quantities appearing in (2) are first evaluated with respect to the centers of the finest boxes and then, if required within the hierarchical algorithm, aggregated to coarser levels. To avoid an expansion of the spatial sources according to (2), we may work directly with the spectral quantitiesJ +M in (3) and represent them by appropriate field expansions within every finest-level box containing equivalent sources.
First, let us consider a distributed SHC expansion according to [1], [19] where the Y nm are the normalized spherical harmonics [20] of degree n and order m according to with the associated Legendre polynomials P m n of degree n and order m. The Y nm is related to the fully 3-D solutions of the scalar Helmholtz equation via [21], [22] Φ nm (r) = z n (kr )Y nm k (6) where z n is an appropriate spherical cylindrical function.
The Cartesian vector f nm contains the unknown expansion coefficients and the dyadic operator T (k) = I −kk with the unit dyadic I converts the Cartesian components into spherical ϑand ϕ-components. Formally, we expand the Cartesian FF components of the magnetic vector potential related to electric volume current densities and convert them into FF electric field components.
As a second possibility, we consider a distributed SHV expansion in the form of where the (spherically transverse) vector functions n nm and m nm are defined as [11], [21], [23] with the spherical unit vectorsê ϑ andê ϕ in ϑand ϕ-direction, respectively. The transverse spherical harmonics are related to the fully 3-D SHVs via The expansion in (7) directly represents the electric FF by a complete and orthonormal vector set of solutions of Maxwell's equations, where c TM nm and c TE nm are the corresponding expansion coefficients.
The lowest order SHCs and SHVs are illustrated in Fig. 2(a) and (b), respectively. Interesting to note is that the sum in (4) starts with n = 0, whereas (7) starts with n = 1. This makes, however, sense, since the dyadic operator T found in (4) partially increases the order of the expansion so that the three dipole fields are obtained for n = 0 in the case of SHC and for n = 1 in the case of SHVs. By the same token, the SHC series is chosen to stop at N − 1 and the SHV series at N. An illustration of a typical distributed expansion is given in Fig. 3. The figure shows a horn alike AUT together with the regular grid of boxes on the finest octree level. Radiation field expansions are assumed in all boxes indicated by a dot, which contain parts of the Huygens' surface around the AUT. The desired Huygens' surface can, e.g., be defined by simple geometric hulls like boxes or ellipsoids, or a mesh can be generated, which helps to identify the non-empty source boxes.

III. HUYGENS' SURFACE ADAPTED DIRECTIVE SPHERICAL HARMONICS EXPANSIONS
The radiation integral contained in (1) is a form of the surface equivalence principle [12, p. 107ff], which states that any field radiated by an AUT outside V AUT can be represented by an equivalent, tangential electric and magnetic surface current density distribution residing on a surrounding Huygens' surface A. If the surface currents on the Huygens' surface are determined in a way that they produce zero fields inside V AUT , then they are directly linked to the tangential fields just outside V AUT according to and we call such surface currents Love currents [14]. Without such a zero-field condition, the surface current densities are not uniquely defined, which is in general not a problem in a regularized iterative or pseudo inverse solution, but removing this kind of redundancy can reduce the number of DOFs in the expansion by half. Imposing the Love condition numerically during the solution of the inverse source problem is possible, but numerically very demanding [5], [13], [18].
As an alternative, the redundancy can also be removed by any other fixed linear relationship between the electric and magnetic surface current densities, where a particularly good choice is an approximation of the Love condition in the form of a combined source condition on the surface A according to [24]  which leads to the concept of the so-called Huygens' radiators [25], [26], [27]. In this way, the number of source DOFs is reduced and the enforcement of such a local condition is computationally very inexpensive. Such combined surface current densities enforce the radiation along the direction of the surface normal vectorn and reduce it toward the opposite direction, i.e., into the AUT volume, where a true null-field can only be achieved for very special configurations. For closed Huygens' surfaces of finite extent, we observe typically a certain field suppression in the exclusion volume, which is, however, rather irrelevant, since our primary aim is just to reduce the number of DOFs while maintaining a good conditioning. Based on these considerations, our goal is to generate DSH expansions starting from (4) and (7), which correspond to equivalent surface current sources according to the combined source condition in (11). The construction of these directive DSHs is performed for the half-space z > 0, i.e., with enforced radiation into this half-space, where it is assumed that equivalent Huygens' surface current densities are only located in the plane z = 0 and the source boxes on the finest level are placed symmetrically around z = 0. Later, the radiation patterns of the obtained basis functions, which are assumed in the non-empty boxes of an appropriate AUT model, as, e.g., illustrated in Fig. 3, are rotated in a way that the former +z-axis points into the direction of the indicated average normals of the boxes. In an efficient implementation of a propagating plane-wave-based inverse source solver, the rotations along ϕ can easily be realized as simple index calculations, if equidistant sampling is assumed in ϕ and if only discrete rotation angles according to integer multiples of the step size in ϕ are assumed. For the rotations in ϑdirection, where a Gauss-Legendre grid [7], [10] is in general more appropriate, the k-space representations of a discrete set of rotated basis functions can be precomputed in the setup phase of the solver.

A. Huygens' Surface Conforming SHVs
We start with an electric surface current density J A defined in the plane z = 0 and verify which modes it can excite in the expansion according to (7). For the c TM/TE nm -coefficients, we find via reciprocity [21] or the bilinear spherical vector mode expansion of the corresponding dyadic Green's functions [22] to be evaluated over the support of J A in the plane z = 0, where the expressions in (8) and (9) have been utilized, and c 1 and c 2 are constants. With the properties of the P nm and its derivatives with respect to ϑ , which need to be evaluated for cos ϑ = 0, i.e., in the plane z = 0, it is found that c TM i.e., both are complementary to each other. From this result, it is immediately obvious that an expansion according to (7), where only the n nm -functions according to (14) and the m nm -functions according to (15) are present, is sufficient to model the fields, since it corresponds to a Huygens' surface integral with electric surface current densities only. Similarly, the dual expansion according to magnetic surface current densities could be used, where the roles of the n nm -and m nm -functions are flipped. However, our goal is to improve the expansion even further and construct an expansion with directive basis functions.
To understand this construction process, let us consider the degree n = 1 first. The n nm -functions of order m = ±1 correspond to circularly polarized Hertzian dipoles with currents oriented in x-and y-directions, whereas the order m = 0 corresponds to a z-oriented Hertzian dipole with linear polarization [28, p. 437]. We first generate Huygens' radiators by combining the TM-and TE-modes of degree n = 1 and m = ±1 according to corresponding to two orthonormal Huygens' radiators, which radiate predominantly into the half-space z > 0. The orders m = −1 and m = +1 indicate left-hand and right-hand circular polarization, respectively. This combination scheme makes obviously sense for all pairs of n nm -and m nm -functions with |m| > 0, but not for m = 0. Based on the duality principle [12, p. 120ff], one can think of the n 10 and m 10 functions as the FF patterns of magnetic and electric loop currents, respectively, with infinitely small radius residing at the origin in the x y-plane. The generation of directive basis functions can, thus, be achieved by combining the FF pattern functions of these loop currents with the FF pattern functions of the corresponding star currents of the other type. An electric star current density can be written as where J L = J Lêϕ refers to the corresponding electric loop current density. The electric FF pattern function of such an electric star current density J S (r ) = I δ(r − R)δ(z)ê ρ of radius R and carrying the current I can be calculated for R → 0 as shown in Appendix A resulting into i.e., its electric field corresponds to the n 20 -function. To achieve a directive FF pattern function, we now combine the n 10 -function corresponding to a magnetic loop current with the n 20 -function corresponding to an electric star current. Similarly, we can combine the fields of an electric loop current with the fields of a magnetic star current, resulting in a summation of the m 10 -and m 20 -functions, respectively. It can be verified that this kind of construction scheme is useful for all n nm -and m nm -functions with |m| < n, where every function just needs to be combined with the corresponding function with degree n + 1 or n − 1. In particular, this combination scheme is mandatory for all m n0 -and n n0 -functions. Based on these observations and by considering the excitation properties in (14) and (15) of our originally assumed electric surface current density in the x y-plane, we propose to construct a first set of directive vector basis functions (SHV-D1) according to where functions including n + 1-terms are only considered if n + 1 ≤ N. With this directive shaping scheme, every excited vector mode according to (14) and (15) is only combined once with another mode, which is not directly excited by an electric surface current density within the x y-plane. If both directive shaping schemes are combined for n > |m| > 0, an alternative second set of directive vector basis functions (SHV-D2) can be constructed according to D (2) nm k = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ nm nm k , |m| = n nm nm k + nm n+1,m k , n > |m| > 0 m n0 k + m n+1,0 k , n odd, m = 0 n n0 k + n n+1,0 k , n even, m = 0 (21)  Fig. 4 up to n = 3. The improved directivity of the SHV-D2 functions as compared to the corresponding SHV-D1 functions is clearly seen. The directivities of the different lowest order SHVs with respect to an isotropic radiator have also been computed and are compared among each other in Table I. The total plane-wave spectrum of every box of the octree containing equivalent sources for the representation of the AUT is finally expanded as where c nm are complex expansion coefficients and where either the functions in (19) or in (21) are used.

B. Huygens' Surface Conforming SHCs-D
Starting with an electric surface current density J A in the plane z = 0, the z-component in the expansion in (4) can immediately be dropped, resulting into a reduction of the DOFs by one third. Next, we verify again which spherical harmonics can be excited by an electric surface current density confined to the plane z = 0. By invoking reciprocity [21] or by working with a bilinear expansion of the scalar Green's function of the Helmholtz equation [22], we find for the x/y-components of the excitation coefficients with Φ nm according to (6). Similar to (14), again with the properties of the P nm to be evaluated for cos ϑ = 0, it is found that only spherical harmonics with n = 0, . . . , N − 1 and m = −n, −(n − 2), . . . , (n − 2), n can be excited, i.e., every other coefficient is zero. Next, Huygens' radiator like expansion functions are created by combining every nonzero electric current mode with a corresponding magnetic current mode. This means nothing where f nm,t ×ê z in the second term forms a magnetic current vector perpendicular to the corresponding electric current vector in the first term, and the operatorê r ×T (k) produces the corresponding electric FF. The subscript t indicates transverse vectors with x/y-components only and the m+=2 in the sum over m indicates that the step size is 2 according to (24). In a more compact notation, the expansion (named SHC-D1) can be written as where the Cartesian-to-spherical transformation operator G is defined as G k = cos(ϑ) cos(ϕ) + cos(ϕ) cos(ϑ) sin(ϕ) + sin(ϕ) cos(ϑ) sin(ϕ) + sin(ϕ) cos(ϑ) cos(ϕ) + cos(ϕ) .
As might be expected, the second directive shaping scheme as discussed for the SHVs with |m| < n can also be applied to the SHCs and, thus, a modified alternative expansion (named SHC-D2) with even better directive shaping behavior can be constructed, where the corresponding spherical harmonics in (26) is replaced according to where again the added functions do not belong to the originally excited set of modes. Some of the SHC-D1 and SHC-D2 basis functions are depicted in Fig. 5, where the improved directive shaping of the SHC-D2 functions as compared to the corresponding SHC-D1 functions is seen for |m| < n. Again, also the directivities of the different SHCs have been computed and are summarized in Table II. The (1, 1)-and (2, 2)-functions even lose some overall directivity due to the directive shaping toward the half-space z > 0, which is due to their reduced directivity within the transverse x y-plane.

IV. NUMERICAL RESULTS
The proposed directive DOF-reduced DSH expansions are evaluated and compared to the corresponding full DSH and equivalent surface current expansions for a variety of NF transformation setups with synthetic and real NF measurement data. Throughout, the generalized minimal residual (GMRES) solver [29] is used for the iterative solution of the resulting normal error systems of equations. The FF obtained from an inverse source solution is compared to another FF via the linear electric field deviation FF (ϑ, ϕ) = 20 log 10 where E FF,1/2 are the FF magnitudes of the two FFs and where the maximum in the denominator is chosen in a certain angular range of interest. If not stated otherwise, a relative stopping criterion rel = 0.997 is used for solver termination, i.e., the solver stops if the residual is reduced by less than rel for two iterations in a row. All computations have been carried out on an Intel 1 Core 2 i7-10700 CPU @ 2.9 GHz workstation with eight cores.

A. Corner Reflector
The first example is a synthetic scattering problem with an open planar Huygens' surface in front of a corner reflector and planar or quasi-planar observations, which have been obtained from a method of moments (MoM) simulation [15], [30] for 1 Registered trademark. 2 Trademarked. Fig. 6. Incident plane-wave normal to the aperture of the corner reflector model. The figure is taken from [24]. The gray-shaded triangle over the aperture is used as Huygens' surface in the inverse source solutions. a plane wave with a frequency of f = 64 GHz normally incident with respect to the aperture. The FF obtained from an inverse source solution corresponds, thus, to the bistatic radar cross section (B-RCS) of the corner reflector, where reference results are available from the MoM solution. One set of NF data is sampled regularly with λ/2 sampling distance in the x y-plane at z = 1.0 m, resulting in 3 243 601 samples for each polarization. A second set of NF data was obtained with irregular sample position distribution, where uniformly distributed random variations in all three spatial directions with respect to the regular samples as used before and with an overall standard deviation of 0.707λ were superimposed [1]. The configuration is depicted in Fig. 6 and due to the planar Huygens' surface it is in particular suitable to investigate the directive behavior of the various equivalent source expansions. The Huygens' surface above the reflector as seen in gray in Fig. 6 is chosen as a triangular surface at z = 0.6 m, which is about 0.02 m above the aperture. For the case of working with surface current densities, the Huygens' surface was meshed with 2 172 301 triangles and RWG functions [8] were used for the inverse source solutions with unconstrained electric and magnetic surface current densities (J + M) or with the related directive Huygens' radiator basis functions (JH) as found in [15]. For DSH expansions, the non-empty source boxes were found as those finest-level octree boxes which are cut by the Huygens' surface, 304 088 in total.
A most relevant B-RCS cut with the main peak due to the triple reflection at ϑ = 0 • and a second peak due to a double reflection at ϑ = −70 • is shown in Fig. 7 for the SHC and SHV inverse source solutions and in Fig. 8 for the corresponding reduced and directive inverse source solutions, both in comparison to the MoM reference results. The valid angular range due to the planar observation surface is about |ϑ| < 75 • . While the forward RCS peak of the MoM solution seen at ϑ = −180 • is correct, there is no need for the inverse source solutions to produce any radiation for |ϑ| > 90 • . The SHV solution shows pretty good suppression of this parasitic radiation, but the SHC solution shows a very strong parasitic forward RCS peak. The reduced and directive solutions in Fig. 8 show both a very good suppression of the parasitic radiation, where the D2 curves show better suppression than the D1 curves in both cases. Kind of interesting is that the SHC-D curves show a strong image of the double-reflection peak, while the SHV-D curves do not, which is an indication Fig. 7. B-RCS cut of the corner reflector in Fig. 6 obtained with full SHC and SHV DSH expansions for irregular NF data. The deviations of both inverse source results are given with respect to the MoM reference. that the SHV-D expansions are somewhat better conditioned than the SHC-D expansions.
The iterative inverse source solution behaviors of the different DSH expansions are illustrated in Figs. 9 and 10. All expansions reach about the same stopping residual, where the convergence for the regular observation data is throughout slightly faster than for the irregular data. Interesting to note is that the decay of the residuals of the reduced directive expansions is in particular faster for the irregular data, where the reduced number of source DOFs is clearly of benefit. Throughout the decay of the residual of the D2 solutions is slightly faster than that of the corresponding D1 solutions. Finally, the solver performances for the different expansions are summarized in Table III. The DSH expansions require throughout less computational resources than the equivalent surface current expansions. The benefits of the reduced directive DSH expansions become especially obvious for the irregular NF data where in particular the number of required iterations is reduced. For recommended "normal-accuracy" settings, the required computation times are about one half of the times given in the table, which are for settings with improved accuracy.

B. DRH18 Double-Ridged Waveguide Antenna
To investigate the introduced DSH expansions for real NF antenna measurements, the DRH18 double-ridged waveguide antenna from RFspin [31] as illustrated in Fig. 11 is considered. The antenna has been measured in the anechoic chamber at the Technical University of Munich for a frequency of 18 GHz. The spherical NF measurements have been performed for two orthogonal polarizations with an open-ended rectangular hollow waveguide probe in a distance of 2.68 m   surface as seen in Fig. 11(b) was meshed with 81 390 triangles, where RWG functions [8] were used for the inverse source solutions with unconstrained electric and magnetic surface current densities (J + M) or with the related directive Huygens' radiator basis functions (JH) as found in [15]. For DSH expansions, the non-empty source boxes were again found as those finest-level octree boxes, which are cut by the Huygens' surface, 2592 in total. A second set of irregular NF observation data was generated out of an inverse source solution of the regular measurement data. A spherical mesh somewhat larger than the minimum sphere of the antenna was here used to avoid too much spatial filtering of the room echoes present in the measurement data. With the obtained sources dual-polarized measurement samples with the same open-ended hollow waveguide probe at 50 000 irregular spherical sample locations with an average distance of 1.5 m from the rotation center were produced. First, a close to uniform distribution of the sample locations was produced according to [32]. Second, random uniformly distributed spatial deviations in all three dimensions with an overall standard deviation of 0.707λ have been added to the regular positions.
The iterative inverse source solution behavior of the SHV and SHC DSH expansions in comparison to the corresponding results of surface current density expansions (J + M and JH) is illustrated in Figs. 12 and 13, respectively. The solution behavior is here mostly determined by the spatial filtering behavior of the expansions, due the relatively large amount of parasitic room echoes within the observation data. The surface current expansions possess of course the best filtering behavior and remove a larger amount of the room echoes, resulting in somewhat larger stopping residuals. For the DSH expansions, the slightly larger residuals of the reduced directive D1 and D2 expansions indicate a slightly better spatial localization behavior than the corresponding SHV and SHC expansions. The irregular NF data leads to slightly larger residuals, probably due to a somewhat improved spatial filtering behavior as compared to the regular data, where the decay rate of the residuals is more or less the same for regular and irregular data. This can be expected, since the influence of the room echoes is here larger than the conditioning effects of the spatial sample distribution.
One FF cut of the antenna is shown in Fig. 14, where reduced directive SHV-D2 solutions for the regular (D2) and for the irregular (D2a) NF data are compared to a solution with unconstrained electric and magnetic surface current densities (J + M). The deviation between the different results is well below the achievable measurement accuracy of the used measurement chamber. NF magnitude distributions in a cut plane through the antenna model as shown in Fig. 11(b) are   surface the suppression effect on the field magnitudes due to the outward oriented directive behavior of the D1 and D2 DSH expansion can be clearly seen, where the D2 expansion show a slightly stronger suppression. Due to the relatively complicated shape of the Huygens' surface and the related arrangement of nonempty octree boxes as, e.g., illustrated in Fig. 3 and an even stronger suppression can certainly not be expected.
An overview of the solver performances of the different expansions is given in Table IV. Again, the DSH Fig. 16. Electric field magnitude obtained with different inverse source solutions in cut a plane through the DRH18 antenna model in Fig. 11(b). expansions need less computational resources than the equivalent surface current expansions. The effect of the irregular data is here, however, less pronounced than in the case of the synthetic data. Consequently, the reduced number of source DOFs of the reduced DSH expansions is here also less pronounced. Due to the rather large room echoes in the data, the spatial filtering properties of the expansions are most relevant for the solution behavior, leading to relatively equal solver performance of all DSH expansions.

V. CONCLUSION
Various DSH expansions for electromagnetic inverse source solutions have been presented and investigated. The expansion centers are placed in the centers of non-empty source boxes of an octree, which is the basis of a hierarchical propagating plane-wave representation of the radiation operator. The nonempty source boxes are identified as those, which are crossed by a suitable Huygens' surface surrounding the relevant source region, where, however, a specific mesh of the Huygens' surface is not necessarily needed. The known standard SHC and SHV DSH expansions are redundant in such an arrangement. Therefore, reduced and directive DSH expansions have been derived by eliminating unnecessary modal functions and by combining different functions to predominantly radiate into the direction of the surface normal on the assumed Huygens' surface away from the AUT volume. Based on the properties of the spherical harmonics, two directive combination schemes have been identified and implemented; the first combines different function types of the same degree and order, the second combines functions of the same type, but with different degrees, where both combination schemes can even be combined. The resulting DSH expansions, two SHV and two SHC expansions based on the different directive combination schemes, reduce the numbers of source DOFs considerably and improve the conditioning of the resulting systems of equations. The performance of the presented DSH expansions has been demonstrated by inverse source solutions for synthetic as well as measured NF observations, where the "doubly-directive" expansions show the best performance throughout.

APPENDIX A FF OF A STAR CURRENT
The FF of a star current J S (r ) = I δ(ρ − R)δ(z )ê ρ of radius R carrying the current I can be calculated via the FF magnetic vector potential A S (r) = μ 0 e −jkr 4πr V J S r e jk r·r r dv (30) where μ 0 denotes the permeability of free space. The dot product r · r with r expressed in spherical coordinates and r in cylindrical coordinates is r · r = ⎛ ⎝ r cos ϕ sin ϑ r sin ϕ sin ϑ r cos ϑ ⎞ ⎠ · ⎛ ⎝ ρ cos ϕ ρ sin ϕ z ⎞ ⎠ = rρ cos ϕ − ϕ sin ϑ + r z cos ϑ.
Inserting this into (30) together with the star current density and integrating over ρ and z with consideration of the filter property of the Dirac delta yields A S (r) = μ 0 I e −jkr 4πr 2π ϕ =0ê ρ e jk R cos(ϕ−ϕ ) sin ϑ R dϕ . (32) Now, we express the unit vectorsê ρ within the integral by Cartesian unit vector aŝ e ρ = cos ϕ ê x + sin ϕ ê y (33) and the Cartesian unit vectors by unprimed cylindrical unit vectors, which are independent of the primed integration variable, resulting intô With this, the integral in (32) can be evaluated as A S (r) = μ 0 I R e −jkr 4πr j2πJ 1 (k R sin ϑ)ê ρ − 0ê ϕ where J 1 is the Bessel function of the first kind. We now consider R → 0 and express the Bessel function accordingly by J 1 (k R sin ϑ) ≈ (k R sin ϑ)/2 resulting into A S (r) = jμ 0 kπR 2 I e −jkr 4πr sin ϑê ρ where the moment πR 2 I of the star current density is assumed to remain constant for R → 0. The FF electric field in spherical coordinates withê ρ = sin ϑê r + cos ϑê ϑ can now be calculated as E S (r) ≈ −jω A S (r) − A S,r (r)ê r = ωμ 0 kπR 2 I e −jkr 4πr sin ϑ cos ϑê ϑ ∼ n 20 . (37)