## 2 Methods
### Numerical Model  
In order to numerically solve the scattering problem in the solution domain \(\mathrm{V}_{int}\), Equations (3) and (4) are cast in weak form as:  
\[\int\limits_{\Sigma}\varphi\Lambda(E^{s})dl-\int\limits_{\mathrm{V}_{int}} \left[\
abla\varphi\cdot\
abla E^{s}-k^{2}m^{2}\varphi E^{s}\right]dS=f, \tag{14}\]  
where \(f\) is a known term depending on the external impinging field:  
\[f=-k^{2}\int\limits_{\mathrm{V}_{int}}\left(m^{2}-1\right)\varphi E^{i}dS, \tag{15}\]  
and Equations (14) and (15) hold for any function \(\varphi\in H^{1}(\mathrm{V}_{int})\), for \(E^{s}\in H^{1}(\mathrm{V}_{int})\), where \(H^{1}(\mathrm{V}_{int})\) is the Sobolev space defined as \(H^{1}(\mathrm{V}_{int})=\left\{\varphi\in L^{2}(\mathrm{V}_{int}):\ \
abla \varphi\in L^{2}(\mathrm{V}_{int})\right\}\), [29]. To obtain the discrete model, first we express the unknown as \(E^{s}(x,y)=\sum_{i=1}^{N}c_{i}\varphi_{i}(x,y)\), where the shape functions are first-order isoparametric nodal elements, the \(c_{i}\)s are the unknown coefficients (degrees of freedom, DoF) and \(N\) is the total number of nodes of the finite element mesh. Then, the Galerkin approach applied to Equation (14) gives the following linear system of equations for the DoFs:  
\[(\mathbf{\Lambda}+\mathbf{A})c=f, \tag{16}\]  
where  
\[\Lambda_{lk}=\int_{\Sigma}\varphi_{l}\Lambda(\varphi_{k})dl,\ \Lambda_{lk}=-\int_{ \mathrm{V}_{int}}\left[\
abla\varphi_{l^{\prime}}\
abla\varphi_{k}-k^{2}m^{2} \varphi_{l}\varphi_{k}\right]dS, \tag{17}\]\[f_{l}=-k^{2}\int\limits_{V_{int}}\!\!\left(m^{2}-1\right)\!\rho_{l}E^{i}dS. \tag{18}\]  
It is convenient partitioning the unknowns DoFs in \(c_{b}\) and \(c_{i}\) which are related to the boundary and the internal nodes of the finite element mesh. In this way, we get:  
\[\left(\left[ \begin{array}{cc}\mathbf{\Lambda}_{bb}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}\end{array} \right]+\left[ \begin{array}{cc}\mathbf{A}_{ bb}&\mathbf{A}_{bi}\\ \mathbf{A}_{ib}&\mathbf{A}_{ii}\end{array} \right]\right)\left[ \begin{array}{c} \mathbf{c}_{b}\\ \mathbf{c}_{i}\end{array} \right]=\left[ \begin{array}{c}\mathbf{f}_{b}\\ \mathbf{f}_{i}\end{array} \right]. \tag{19}\]  
We notice that \(\mathbf{\Lambda}_{bb}\) is a \(N_{b}\times N_{b}\) fully populated matrix, whereas \(\mathbf{A}\) is a \(N\times N\) sparse matrix, where \(N_{b}\) is the number of boundary nodes, usually much smaller than \(N\) (\(N_{i}\)). Matrix \(A\) is the classical stiffness matrix arising from FEM. For solving the linear system (19), an iterative method is preferable to a direct method. In doing this, there are several key aspects such as memory occupation and computational cost for a matrix-by-vector product. The memory required to store \(\mathbf{\Lambda}\) and \(\mathbf{A}\) is \(N_{b}^{2}\) and \(O(N)\), respectively. On the other hand, the number of multiplications for a matrix-by-vector product is \(N_{b}^{2}+O(N)\), which can be further reduced by properly treating the stiffness matrix. To make a comparison with the PML, we notice that the memory occupation and the number of multiplications is \(O(N_{b})+O(N_{PML})\), where \(N_{PML}\) are the number of DoFs in the computational domain for the PML and \(N_{b}\) the number of DOFs located in the absorbing layer of thickness \(\delta\), see Equation (8). As mentioned, the computational domain for the PML has to be far from the scatterer, thus it should be \(N_{PML}\gg N\), for a prescribed mesh density. Regarding iterative methods, in this work we applied the biconjugate gradient technique [30], with an Incomplete Lower-Upper (ILU) factorization used as a preconditioner [31]. In conclusion, at a given accuracy, the DtN approach, if compared to the PML one, does not introduce extra DoFs in the PML layer (\(N_{b}\)) and can strongly reduce the inner DoFs (\(N_{i}\)), since the boundary of the computational domain can be placed close to the scatterers.