## 2 Formulation of the Problem and Proposed Solution
Taking advantage of the revolution symmetry of the problem at hand, the fields can be expanded in Fourier series, i.e.,  
\[f(\rho,\phi,z)=\sum_{n=-\infty}^{+\infty}f^{(n)}(\rho,z)e^{jnp}. \tag{2}\]  
In this way, the orthogonality of the Fourier series terms allows to reduce the system of surface integral equations to an infinite set of independent systems of 1-D singular integral equations for the harmonics of the effective surface current densities. It is possible to show that the system of equations associated to the \(n\)-th harmonics of the currents can be written in the VHT domain as in the following:  
\[\sum_{j=1}^{2}\int\limits_{0}^{+\infty}\underline{\mathbf{H}}^{(n)}(wp) \bigg{(}\underline{\widetilde{\mathbf{G}}}(w)-\delta_{i,j}R_{\underline{ \mathbf{I}}}\bigg{)}\underline{\widetilde{\mathbf{J}}}_{j}^{(n)}(w)e^{- \mathrm{i}|z_{i}-z_{j}|\sqrt{k_{0}^{2}-w^{2}}}wdw=-\underline{\mathbf{E}}^{ \mathrm{inc}(n)}(\rho,z_{i}) \tag{3}\]  
for \(\rho\leq a\) and \(i=1,2\), where the symbol \(\underline{\mathbf{P}}^{(n)}(\cdot)\) is introduced to denote the column vector \(\left( \begin{array}{c}P_{\rho}^{(n)}(\cdot)\\ -\mathrm{i}P_{\phi}^{(n)}(\cdot)\end{array} \right)\),  
\[\underline{\widetilde{\mathbf{J}}}_{i}^{(n)}(w)=\mathrm{VHT}_{n}\left[ \underline{\widetilde{\mathbf{J}}}_{i}^{(n)}(\rho)\right]=\int\limits_{0}^{+ \infty}\underline{\mathbf{H}}^{(n)}(w\rho)\underline{\widetilde{\mathbf{J}}} _{i}^{(n)}(\rho)\rho d\rho \tag{4}\]  
is the VHT of order \(n\) (VHT\({}_{n}\)) of the \(n\)-th harmonic of the current on the \(i\)-th disk, \(\underline{\mathbf{H}}^{(n)}(w\rho)\) is the kernel of the VHT\({}_{n}\)[44],  
\[\underline{\widetilde{\mathbf{G}}}(w)=\left( \begin{array}{cc}\widetilde{G}_ {CC}(w)&0\\ 0&\widetilde{G}_{DD}(w)\end{array} \right)=\frac{1}{2\omega\varepsilon_{0}} \left( \begin{array}{cc}-\sqrt{k_{0}^{2}-w^{2}}&0\\ 0&-k_{0}^{2}/\sqrt{k_{0}^{2}-w^{2}}\end{array} \right) \tag{5}\]  
with \(\sqrt{k_{0}^{2}-w^{2}}=-\mathrm{i}\,\sqrt{-k_{0}^{2}+w^{2}}\) is the spectral domain Green's function of the problem at hand, \(\delta_{i,j}\) is the Kronecker delta, \(\underline{\mathbf{I}}\) is the identity matrix, and \(\underline{\mathbf{E}}^{\mathrm{inc}(n)}(\rho,z_{i})\) is the \(n\)-th harmonic of the incident electric field in the \(i\)-th disk plane. By means of the Helmholtz decomposition, the \(n\)-th harmonic of the current on the \(i\)-th disk is replaced by the superposition of a surface curl-free contribution, \(\underline{\widetilde{\mathbf{J}}}_{i,C}^{(n)}(\rho)\), and a surface divergence-free contribution, \(\underline{\widetilde{\mathbf{J}}}_{i,D}^{(n)}(\rho)\)[45]. This is an advantageous choice because it allows us to obtain scalar unknowns in the spectral domain. Indeed, it is simple to show that the VHT\({}_{n}\) of \(\underline{\widetilde{\mathbf{J}}}_{i,T}^{(n)}(\rho)\) with \(T=\mathrm{C},\mathrm{D}\) are scalar functions [41], i.e.,  
\[\mathrm{VHT}_{n}\left[\underline{\widetilde{\mathbf{J}}}_{i,C}^{(n)}(\rho) \right]=\left( \begin{array}{c}\widetilde{\widetilde{\mathbf{J}}}_{i,C}^{(n) }(w)\\ 0\end{array} \right), \tag{6}\]  
\[\mathrm{VHT}_{n}\left[\underline{\widetilde{\mathbf{J}}}_{i,D}^{(n)}(\rho) \right]=\left( \begin{array}{c}0\\ -\mathrm{j}\widetilde{\mathbf{J}}_{i,D}^{(n)}(w)\end{array} \right). \tag{7}\]  
Now, the system of Equation (3), which cannot admit a closed form solution, has to be solved by resorting to a discretization scheme. It is important to note that fast convergence and analytical regularization can be simultaneously achieved by means of the Galerkin method providing suitable expansions of the unknowns. In order to guarantee fast convergence, orthogonal bases of the functional spaces at which the unknowns belong are the best choice. Such spaces can be completely characterized by the edge behavior [47] and the behavior around the origin of the harmonics of the currents, and observing that the sources are off the disks:  
\[J_{i,t}^{(n)}(\rho)= \begin{cases}\rho^{||n|-1}|(a-\rho)^{p_{\text{T}}}\overline {J}_{i,t}^{(n)}(\rho)&\rho<a\\ 0&\rho>a\end{cases}