## 3 Fast Evaluation of the Matrix Coefficients  
For \(n,m\geq 0\), the general integral associated with the primary contribution can be written as  
\[I_{P}=\int\limits_{C}f_{i}\big{(}u,\overline{x}_{j},\overline{z}_{j}\big{)} \widetilde{\phi}_{n}^{(\alpha_{i},\beta_{i})}(a_{i}u)\widetilde{\phi}_{m}^{( \alpha_{i},\beta_{j})}\big{(}-a_{j}\big{[}T_{ij}\big{(}u,\overline{x}_{j}, \overline{z}_{j}\big{)}\big{]}^{*}\big{)}du \tag{7}\]  
where \(f_{i}(u,\overline{x}_{j},\overline{z}_{j})\) is one of the kernels in Equations (4a)-(4d), \(T_{ij}(u,\overline{x}_{j},\overline{z}_{j})\) is defined in Equation (A1f), and the star denotes the complex conjugate of a complex number, while the general integral associated to the scattered contribution is  
\[I_{S}=I_{S}^{\prime}+I_{S}^{\prime\prime} \tag{8a}\] \[I_{S}^{\prime}=\int\limits_{C}\big{[}g_{i}\big{(}u,\overline{x}_{j}, \overline{z}_{j}\big{)}-g_{i}^{\infty}\big{(}u,\overline{x}_{j},\overline{z}_ {j}\big{)}\big{]}\widetilde{\phi}_{n}^{(\alpha_{i},\beta_{i})}(a_{i}S_{i}(u)) \widetilde{\phi}_{m}^{(\alpha_{i},\beta_{i})}\big{(}-a_{j}\big{[}S_{j}(u)\big{]} ^{*}\big{)}du\] (8b) \[I_{S}^{\prime\prime}=\int\limits_{C}g_{i}^{\infty}\big{(}u,\overline{x}_{j}, \overline{z}_{j}\big{)}\widetilde{\phi}_{n}^{(\alpha_{i},\beta_{i})}(a_{i}S_{i }(u))\widetilde{\phi}_{m}^{(\alpha_{i},\beta_{j})}\big{(}-a_{j}\big{[}S_{j}(u) \big{]}^{*}\big{)}du \tag{8c}\]  
where \(g_{i}\big{(}u,\overline{x}_{j},\overline{z}_{j}\big{)}\) is one of the kernels in Equations (4e)-(4h), \(g_{i}^{\infty}\big{(}u,\overline{x}_{j},\overline{z}_{j}\big{)}\) is obtained starting from \(g_{i}\big{(}u,\overline{x}_{j},\overline{z}_{j}\big{)}\) by replacing \(C_{K}(u)\) for \(K\in\{E,H\}\), defined in Equations (A1m) and (A1n), respectively, with \(C_{K}^{\infty}(u)=\lim\limits_{d\to+\infty}C_{K}(u)\), and \(S_{j}(u)\) is defined in Equation (A1p). Since the coefficient matrix elements for \(n=-1\) and/or \(m=-1\) can be written in a similar way, they are here omitted for the sake of brevity. The integrals in Equations (7), (8b) and (8c) are improper integrals. According to the asymptotic behavior for the large argument of the CHF1 [45], i.e.,  
\[{}_{1}F_{1}(a;b;u)\overset{|u|+\to+\infty}{\sim}\frac{e^{\pm j\pi a}u^{-a} \Gamma(b)}{\Gamma(b-a)}+\frac{e^{u}u^{a-b}\Gamma(b)}{\Gamma(a)} \tag{9}\]  
where the upper sign has to be chosen for \(-\pi/2\leq\arg(u)<3\pi/2\) and the lower sign for \(-3\pi/2<\arg(u)\leq\pi/2\), and observing that \(C_{K}(u)-C_{K}^{\infty}(u)\overset{|u|+\to+\infty}{\sim}e^{-2|u|d}\), it can be concluded that the integrand of the integral in Equation (8b) has an asymptotic exponential decay, while the integrals in Equations (7) and (8c) involve asymptotically oscillating functions with an algebraic decay in the worst cases. As a result, the numerical evaluation of such kinds of integrals is less and less efficient as the accuracy required for the solution is higher. The asymptotic decay of the integrands in Equations (7) and (8c) can be sped up by pulling out the asymptotic behavior of the integrand such that the integrals of the extracted contributions can be expressed in closed form. Nevertheless, such a procedure does not change the oscillating nature of the integrands, which keep on having an algebraic decay. Therefore, the numerical convergence of the accelerated integrals, even if faster, is still strongly dependent on the accuracy required for the solution. A new procedure, overcoming this problem, is shown in the following. For the sake of brevity, only the integrals associated to the primary contribution for \(n,m\geq 0\) and the propagation of bound modes are considered, the generalization of the proposed technique to the other cases being straightforward. For bound modes propagation, the propagation constant is real, i.e., \(v=v_{R}\), and such that \(k_{2}<v_{R}<k_{1}\); therefore, the square-roots branch points of the free-space Green's function in the spectral domain are imaginary numbers.