## 2 Formulation of the Problem and Proposed Solution
### Proposed Solution
\begin{bmatrix}-\mathbf{j}2\omega\epsilon_{2}\mathbf{A}_{ \mathrm{C},\mathrm{C}}^{(n)}&-\mathbf{j}\sqrt{\frac{\epsilon_{2}}{\mu_{2}}} \mathbf{A}_{\mathrm{C},\mathrm{D}}^{(n)}\\ 2\mathbf{j}\sqrt{\frac{\epsilon_{2}}{\mu_{2}}}\mathbf{A}_{\mathrm{D},\mathrm{ C}}^{(n)}&\mathbf{j}\frac{2}{\omega\epsilon_{2}}\mathbf{A}_{\mathrm{D}, \mathrm{D}}^{(n)}\end{bmatrix} , \tag{20a}\] \[\mathbf{A}_{\mathrm{C},\mathrm{C}}^{(n)}=\left\{A_{\mathrm{C},\mathrm{C},k,h}^{(n) }\right\}_{k=-1+\delta_{n,0}}^{+\infty},\] (20b) \[\mathbf{A}_{\mathrm{C},\mathrm{D}}^{(n)}=\left\{A_{\mathrm{C},\mathrm{D},k,h}^{(n )}\right\}_{k=-1+\delta_{n,0},\ h=0}^{+\infty},\] (20c) \[\mathbf{A}_{\mathrm{D},\mathrm{C}}^{(n)}=\left\{A_{\mathrm{D},\mathrm{C},k,h}^{(n )}\right\}_{k=0,\ h=-1+\delta_{n,0}}^{+\infty},\] (20d) \[\mathbf{A}_{\mathrm{D},\mathrm{D}}^{(n)}=\left\{A_{\mathrm{D},\mathrm{D},k,h}^{(n )}\right\}_{k=0}^{+\infty},\] (20e) \[A_{\mathrm{C},\mathrm{C},-1,-1}^{(n)}=\left(\hat{\alpha}^{(n)}\right)^{2}\left[M _{\mathrm{C},-1,-1}^{(n)}-\left(\alpha^{(n)}\right)^{2}M_{\mathrm{D},-1,-1}^ {(n)}\right]\ \ \mathrm{for}\ \ \ n\
eq 0,\] (20f) \[A_{\mathrm{C},\mathrm{C},-1,h}^{(n)}=A_{\mathrm{C},\mathrm{C},h,-1}^{(n)}=\hat{ \alpha}^{(n)}M_{\mathrm{C},-1,h}^{(n)}\ \ \mathrm{for}\ \ \ n\
eq 0,\ h\geq 0,\] (20g) \[A_{\mathrm{C},\mathrm{C},k,h}^{(n)}=A_{\mathrm{C},\mathrm{C},k,h}^{(n)}=M_{ \mathrm{C},k,h}^{(n)}\ \ \mathrm{for}\ \ k,h\geq 0, \tag{20h}\]  
integrals in (21f) and (21g) are asymptotically exponentially decaying functions, the corresponding integrals quickly converge. On the other hand, the first integral in (20l) and the first integral in (20m), which are improper integrals of oscillating functions with an asymptotic algebraic decay, can be rewritten as combinations of fast-converging proper integrals by means of an analytical procedure in the complex plane proposed by the authors in a previous paper [26]. Henceforth, the numerical simulations were provided by means of in-house C++ software code running on a laptop equipped with an Intel Core i7-10510U 1.8 GHz processor and 16 GB of RAM, and the integrals were evaluated by means of an adaptive scheme based on a 24-point Gauss-Legendre quadrature rule. The fast convergence of the proposed method can be quantified by introducing the following normalized truncation error:  
\[\text{err}(M,N)=\sqrt{\frac{\sum\limits_{n=-N+1}^{N-1}\left\|\mathbf{x}_{M+1}^ {(n)}-\mathbf{x}_{M}^{(n)}\right\|^{2}}{\sum\limits_{n=-N+1}^{N-1}\left\| \mathbf{x}_{M}^{(n)}\right\|^{2}}}, \tag{22}\]  
where the symbols \(\left\|\cdot\right\|\) and \(\mathbf{x}_{M}^{(n)}\) denote, respectively, the Euclidean norm and the vector of the first \(M\) expansion coefficients for each unknown, whereas \(2N-1\) is the number of azimuthal harmonics used. It can be shown that the convergence rate is slightly affected by the disk depth and the position and orientation of the dipole. However, as expected, it is affected by the disk radius in terms of both the number of expansion functions and azimuthal harmonics to be used in order to achieve a predetermined accuracy. Such behavior is shown in Figure 2,b, where the normalized truncation error is plotted as a function of \(M\) for \(N=7\) and as a function of \(N\) for \(M=8\) for a zero-thickness PEC disk of radius \(a=5,\;10,\;15,\;20\) cm buried in a lossy half-space with \(\varepsilon_{r2}=3.5-\text{j}0.3\) and \(\mu_{r2}=1\) at a depth \(d=3\) cm, and illuminated by a Hertzian dipole directed along the \(y\) axis with \(\left|I_{0}\right|=10^{-3}\text{Am}\), operating at 800 MHz and located in the upper half-space filled by air, \(\varepsilon_{r1}=\mu_{r1}=1\), at \(x_{0}=-2\) m, \(y_{0}=2\) m, \(z_{0}=2\) m. It is worth noting that \(M=8\) and \(N=7\) achieved a normalized truncation error less than \(10^{-3}\) for the worst case examined, i.e., for \(a=20\) cm, with a computation time of 50 s.  
Validation of the developed in-house C++ software code was carried out by means of comparisons with CST-MWS.