## 3 Proposed Solution  
The Helmholtz-Galerkin technique is adopted to simultaneously regularize and discretize the dual integral Equation (14).  
The surface curl-free and divergence-free contributions of \(\mathbf{j}_{i}^{(n)}(\rho)\) are assumed to be new unknowns in the spatial domain, thus leading to the scalar unknowns \(\tilde{J}_{i,\mathrm{C}}^{(n)}(w)\) and \(\tilde{J}_{i,\mathrm{D}}^{(n)}(w)\) in the spectral domain [36].  
The behavior of \(\mathbf{j}_{i}^{(n)}(\rho)\) at the disk rim and around the center of the disk can be reconstructed by expanding the unknowns in the spectral domain in the following complete and non-redundant series of weighted Bessel functions of the first kind [35; 56]:  
\[\tilde{J}_{i,T}^{(n)}(w)=\sum\limits_{h=-1+\delta_{n,0}}^{+\infty}\gamma_{T,h} ^{(n)}\tilde{J}_{h}^{(n)}(w), \tag{16a}\] \[\tilde{J}_{h}^{(n)}(w)=\sqrt{2(|n|+2h+2)}\frac{J_{|n|+2h+2}(aw)}{w}, \tag{16b}\]  
where \(T=\mathrm{C},\mathrm{D}\), \(\delta_{n,m}\) is the Kronecker delta function and \(\gamma_{T,h}^{(n)}\) are the expansion coefficients, orthonormal over the interval \((0,+\infty)\) with the weight function \(w\)[57], i.e.,  
\[\int\limits_{0}^{+\infty}\tilde{J}_{k}^{(n)}(w)\tilde{J}_{h}^{(n)}(w)wdw=2 \sqrt{(|n|+2k+2)(|n|+2h+2)}\int\limits_{0}^{+\infty}\frac{1}{w}J_{|n|+2k+2}( aw)J_{|n|+2h+2}(aw)dw=\delta_{h,k} \tag{17}\]  
with \(h,k=-1+\delta_{n,0},\ldots\). Moreover, the condition \(\gamma_{\mathrm{D},-1}^{(n)}=\mathrm{j}\,\mathrm{sgn}(n)\gamma_{\mathrm{C}-1} ^{(n)}\) for \(n\
eq 0\), where \(\mathrm{sgn}(\cdot)\) denotes the Signum function [55], has to be imposed in order to guarantee the vanishing of \(\mathbf{j}_{i}^{(n)}(\rho)\) for \(\rho>a\), i.e., in order to directly satisfy the equation in the external domain (14b).  
The following matrix equation is obtained by applying the [PERSON] method:  
\[\sum\limits_{h=-1+\delta_{n,0}}^{+\infty}\gamma_{T,h}^{(n)}\int\limits_{0}^{ +\infty}\frac{\widetilde{G}_{T}(w)-R_{i}}{\widetilde{G}_{T}(w)-R_{e}}\tilde{ J}_{k}^{(n)}(w)\tilde{J}_{h}^{(n)}(w)wdw=-\mathrm{j}^{n+1}e^{-\mathrm{i}\mu \phi_{0}}E_{0t}\frac{\tilde{J}_{k}^{(n)}\big{(}k_{0}s_{\theta_{0}}\big{)}}{ \widetilde{G}_{T}\big{(}k_{0}s_{\theta_{0}}\big{)}-R_{e}} \tag{18}\]  
with \(k=-1+\delta_{n,0},\ldots\). Starting from (15) and (17), it is simple to show that  
\[\int\limits_{0}^{+\infty}\frac{\widetilde{G}_{T}(w)-R_{i}}{\widetilde{G}_{T}( w)-R_{e}}\tilde{J}_{k}^{(n)}(w)f_{h}^{(n)}(w)wdw=\int\limits_{0}^{+\infty} \underbrace{\left(\frac{\widetilde{G}_{T}(w)-R_{i}}{\widetilde{G}_{T}(w)-R_{e }}-\kappa_{T}\right)}_{\mathrm{O}(1/w)}\tilde{J}_{k}^{(n)}(w)f_{h}^{(n)}(w)wdw +\kappa_{T}\delta_{h,k}, \tag{19}\]  
where \(\kappa_{\mathrm{C}}=1\) and \(\kappa_{\mathrm{D}}=R_{i}/R_{e}\). Hence, the most singular part of the integral operator is diagonalized. Moreover, following the line of reasoning in [35], it is possible to demonstrate that the remaining part of the discretized operator is a compact operator in \(l^{2}\). Since the discretized counterpart of the free term has a bounded \(l^{2}\)-norm (the sources are off the scatterer surface), the obtained matrix equation is of the Fredholm second kind in \(l^{2}\), i.e., the convergence is guaranteed. On the other hand, the reconstruction of the physical behavior of the unknowns ensures the fast convergence of the method, i.e., few expansion functions lead to highly accurate results. Additionally, the matrix coefficients, which are improper integrals involving products of Bessel functions, are efficiently numerically evaluated by means of suitable analytical techniques recently developed by the authors [35; 52].