## 2 Formulation of the Problem and Proposed Solution
, \tag{8}\]  
for \(t=\rho,\phi\), \(p_{\rho}=1/2\), \(p_{\phi}=0\), where \(\overline{J}_{i,t}^{(n)}(\rho)\) are well behaved functions. By selecting the following expansion functions:  
\[\widetilde{f}_{T,h}^{(n)}(w)=\sqrt{2(|n|+2h+p_{\text{T}}+1)}\frac{J_{|n|+2h+p_ {\text{T}}+1}(aw)}{w^{p_{\text{T}}}}, \tag{9}\]  
where \(J_{n}(\cdot)\) is the Bessel function of the first kind [48], which are orthonormal functions on the interval \((0,+\infty)\) with the weight function \(w^{2p_{\text{T}}-1}\)[49], it is possible to demonstrate that the complete and non-redundant Neumann series [50]  
\[\widetilde{J}_{i,\text{T}}^{(n)}(w)=\sum\limits_{h=-1+\delta_{k,0}}^{+\infty} \gamma_{i,T,h}^{(n)}\ \widehat{\gamma}_{T,h}^{(n)}(w), \tag{10}\]  
where \(\gamma_{i,T,h}^{(n)}\) is the general expansion coefficient, \(\gamma_{i,\text{D}-1}^{(n)}=\alpha^{(n)}\gamma_{i,\text{C}-1}^{(n)},\alpha^{ (n)}=\frac{j\text{sgn}(n)|n|! }{\Gamma(|n|+3/2)}\sqrt{\frac{a(|n|+1/2)}{2|n|}}\), \(p_{\text{C}}=3/2\) and \(p_{\text{D}}=1\), reconstruct the physical behavior in (8) [40]. On the other hand, by extracting the asymptotic behavior of the kernel for \(i=j\) in (3),  
\[\underline{\widetilde{\mathbf{G}}}(w)-R_{i}\underline{\underline{\mathbf{I}} }\stackrel{{ w\rightarrow+\infty}}{{\sim}}\underline{\widetilde{ \mathbf{G}}}(w)=\left( \begin{array}{cc}\widetilde{\widetilde{\mathbf{G}}}_{ CC,j}(w)&0\\ 0&\widetilde{\widetilde{G}}_{DD,j}(w)\end{array} \right)=\left( \begin{array}{cc}j \frac{w}{2\omega\omega_{0}}&0\\ 0&-R_{i}\end{array} \right), \tag{11}\]  
the system of Equation (3) can be rewritten as in the following:  
\[ \begin{split}&\int\limits_{0}^{+\infty}\underline{\underline{ \mathbf{H}}}^{(n)}(w\rho)\underline{\underline{\mathbf{G}}}_{i}(w)\underline{ \widetilde{\mathbf{I}}}_{i}^{(n)}(w)wdw+\\ &+\sum\limits_{j=1}^{2}\int\limits_{0}^{+\infty}\underline{\underline{ \mathbf{H}}}^{(n)}(w\rho)\left[\underline{\widetilde{\mathbf{G}}}(w)-\delta _{i,j}\left(R_{j}\underline{\underline{\mathbf{I}}}+\underline{\underline{ \mathbf{G}}}_{i}(w)\right)\right]\underline{\widetilde{\mathbf{I}}}_{i}^{(n)}(w )e^{-\|z_{i}-z_{j}|\sqrt{\sqrt{k_{0}^{2}-w^{2}}}}wdw=\\ &=-\underline{\mathbf{E}}^{\text{inc}(n)}(\rho,z_{i})\end{split}  \tag{12}\]  
for \(\rho\leq a\) and \(i=1,2\). The discretization of (12) by means of the Galerkin method with the expansions (10) leads to a matrix operator equation. It is simple to show that the matrix operator associated with the first term at the left-hand side of (12) is diagonal due to the orthonormality property of the expansion functions. The general element of the matrix operator \(\underline{\underline{A}}^{(n)}\) associated with the second term at the left-hand side of (12) can be expressed as a linear combination of integrals of the kind  
\[\int\limits_{0}^{+\infty}\widehat{f}_{T,k}^{(n)}(w)\left[\widetilde{G}_{ \text{TT}}(w)-\delta_{i,j}\left(R_{i}-\widetilde{\widetilde{G}}_{TT,i}(w) \right)\right]\widehat{f}_{T,h}^{(n)}(w)e^{-\|z_{i}-z_{j}|\sqrt{k_{0}^{2}-w^{2} }}wdw. \tag{13}\]  
Looking at Formula (13), it is clear as the integrands of the integrals associated to the mutual contributions, i.e., for \(i\
eq j\), have an asymptotic exponential decay related to the distance between the disks, \(|z_{i}-z_{j}|\). On the other hand, the elements for \(i=j\) coincide with the ones obtained for a single resistive disk analyzed in [40]. Hence, following the same line of reasoning presented in [40], it can be shown that \(\underline{A}^{(n)}\) is a compact operator in \(l^{2}\). To conclude, the obtained matrix equation is of the Fredholm second-kind in \(l^{2}\) because the free term has a bounded \(l^{2}\)-norm [39; 40].