## 1 Introduction  
Computational Electromagnetics techniques have traditionally been utilized in the mathematical modeling of problems related to radiation and the propagation of electromagnetic waves. Standard applications involve antennas, diffraction, scattering, waveguides, propagation in complex environments, etc., while boundary surfaces are typically either perfectly conducting or dielectric. Recently, important technological advances have been based on the electromagnetic properties of Near-Zero Index (NZI) materials [1; 2], plasmonics [3], metasurfaces [4], graphene [5], nanoparticles [6], etc. Analytical and asymptotic methods were originally the only mathematical tool to perform relevant calculations, naturally being restricted to canonical geometries. However, computers demonstrated explosive progress over the last few decades, thus rendering numerical techniques indispensable in extracting approximate results of controllable accuracy for arbitrary geometries. In general terms, numerical techniques discretize either integral [7; 8] or differential equations [9; 10] involved in the mathematical simulation of physical problems. A variant of differential equation techniques is based on the construction of equivalent circuit or transmission line models [11; 12; 13]. Each group of the aforementioned techniques is characterized by advantages and disadvantages, and therefore their suitability and efficiency depend on the particular configuration. For example, in scattering problems, where the Radar Cross Section (RCS) of a target must be calculated, differential equation methods should ideally discretize the entire space all the way to infinity, which is obviously impossible. Mesh truncation is a complicated task that needs to be resolved in order to obtain accurate results. On the contrary, integral equation methods are usually more complicated in terms of calculations, and the computational cost may be exceedingly high due to dense matrices involved in the resulting linear system of equations. However, field behavior at infinity is automatically incorporated into the Green's function, and therefore modeling the entire space is easily tractable. Integral equation solutions and their procedure optimization is an important research topic, which recently triggered some especially unconventional approaches [14]. The method of auxiliary sources (MAS) [15] is an integral equation technique that has successfully been used in computational physics, including a wide range of electromagnetic radiation and scattering applications [16]. MAS is superficially similar to the point matching version of the method of moments (MoM) [7]; however, the auxiliary current sources are located at a distance from, and not right on the surface boundaries. Moreover, the MAS basis functions set, used in the field expansions, has been proven to be complete [17], which is not always easy to prove in MoM, therefore the method is mathematically rigorous. Finally, unlike MoM, MAS is free of singularity complications, it does not require any time-consuming numerical integrations and is much easier to implement algorithmically. A quantitative assessment of MAS performance is given in [18], where a detailed comparison is carried out between the computational costs of MAS and MoM. It is demonstrated that for the same discretization density, MAS is always more efficient than MoM, but even if a higher number of unknowns is required for MAS, the latter is still less computationally intensive than MoM under certain conditions. Although MAS has been invoked in several problems with various configurations and material properties, further research is necessary to determine the optimal source location for arbitrary geometry layouts. Particular difficulties arise when the outer boundary of the scatterer contains wedges, i.e., when the analytical expression of the boundary curve is not differentiable. In that case, the solution accuracy deteriorates because the boundary condition close to the wedge tip is hard to satisfy sufficiently well.