## 2 Methods  
### Scattering Problem and Helmholtz Equation  
The reference problem for this paper is depicted in Figure 1a, where an electromagnetic wave hits a number of objects, giving rise to a scattered field. The scatterers are dielectric but non-magnetic objects. Hereafter we assume that the scatterers are invariant along the \(z\) axis. At any spatial position, \(\mathbf{r}\), and at any time instant, \(t\), the total field is given by the sum of the two aforementioned contributions:  
\[\mathbf{E}(\mathbf{r},t)=\mathbf{E}^{i}(\mathbf{r},t)+\mathbf{E}^{i}(\mathbf{r},t). \tag{1}\]  
The evaluation of the scattering from objects is a classic electromagnetic problem that may be formulated, for instance, by using Helmholtz equations [11]. To this end, in the following we consider time-harmonic fields,  
\[\mathbf{e}(\mathbf{r},t)=\text{Re}\Big{\{}\mathbf{E}(\mathbf{r})e^{i\omega t}\Big{\}}, \tag{2}\]  
where \(\omega\) is the angular frequency. Hereafter, we assume \(\text{TM}_{\text{z}}\) propagation, that is the incident wave is a Transverse Magnetic (TM) field with the electric field polarized along the \(z\) axis and the scatterer is invariant in the \(z\) direction. Therefore, the scattering problems reduces to a 2D scalar problem, where \(E^{i}\) is the (scalar) unknown. The scattered field is the solution of coupled Helmholtz equations written inside and outside the domain \(\Omega\), defined as the union of all the scatterers \(\Omega_{i}\) (\(\Omega=\Omega_{1}\cup\Omega_{2}\cup\ldots\)), see Figure 1a. Proper boundary conditions and regularity conditions at infinity must be imposed to close the problem. The final model is:  
\[\
abla^{2}E+k^{2}m^{2}E=0,\ in\ \Omega, \tag{3}\]  
\[\
abla^{2}E+k^{2}E=0,\ in\ R^{3}\backslash\Omega, \tag{4}\]  
\[E|_{\partial\Omega^{+}}-E|_{\partial\Omega^{-}}=0, \tag{5}\]  
\[\partial_{n}E|_{\partial\Omega^{+}}-\partial_{n}E|_{\partial\Omega^{-}}=0, \tag{6}\]  
Figure 1: The analyzed scattering problem: \(n\) arbitrarily shaped objects are illuminated by a known incident field. (**a**) References for the Helmholtz problem; (**b**) references for the Dirichlet-to-Neumann (DhN) formulation. \[\lim_{r\to\infty}r^{1/2}(\partial_{n}E^{s}-jkE^{s})=0. \tag{7}\]  
Here, \(k\) is the free-space wavenumber, \(m(\mathbf{r})\) is the refraction index of the scatterers, and \(\partial_{n}\) is the outward normal derivative. Specifically, conditions (5) and (6) enforce the continuity of the tangential components of the electrical and magnetic fields at the boundary of the scatterers, \(\partial\Omega\), whereas (7) is the classical Sommerfeld radiation condition. It is worth noting that the problem (Equations (3)-(7)) is defined in an unbounded domain. The numerical solution of the scattering problem (Equations (3)-(7)) may be obtained from integral or differential formulations [10; 11]. The main advantage of using the integral formulations is the reduced solution domain, which includes only the material regions. Conversely, a weak point is the need of the Green function, which is available in analytical form only for simple geometries and material property distributions. In the general case, its numerical evaluation is extremely cumbersome because of the presence of singular or hyper-singular kernels. In addition, integral formulations give rise to a final numerical model characterized by fully populated matrices, which may pose challenging problems for storage and inversion. Conversely, differential formulations do not require knowledge of the Green function and give rise to sparse matrices. However, a weak point is that the solution must be sought in the whole space. To cope with this problem, an artificial boundary, \(\Sigma\), is usually introduced, in order to obtain a bounded solution domain, \(V_{int}\), that includes all the scatterers (Figure 1b). In propagation problems, the presence of this artificial boundary may affect the solution by introducing spurious non-physical reflections.