## 3 Conclusions  
In this work, we considered a fractional Helmholtz equation approximated by ad hoc centered differences with variable wave number \(\mu(x,y)\), in the specific case where the complex-valued function \(\mu(x,y)\) has a pole of order \(\gamma\). Eigenvalue distribution and clustering results have been derived in [1; 2]. The numerical results presented in this corroborate the analysis.  
Many more intricate cases can be treated using the same type of theoretical apparatus, including the GLT theory [3; 18] and non-Hermitian perturbation results, such as those in [4; 6]. We list a few of them.  
* The numerical results in Section 2.5 seem to indicate that the spectral distribution of the original matrix sequence and the spectral clustering at \(1\) of the preconditioned matrix sequence holds also when the Frobenius norm condition in [4] is violated; this is an indication that Theorem 1 in [4] may not be sharp. A related conjecture is that the key condition \(\|Y_{n}\|_{2}^{2}=o(n)\) in Theorem 2 could be replaced by \(\|Y_{n}\|_{p}^{p}=o(n)\), with any \(p\in[1,\infty)\), which would be very useful when the trace norm is considered, i.e., for \(p=1\).
* Definition 1 has been reported with a matrix size of the symbol equal to \(r\geq 1\). In our study for matrices arising from finite differences, the parameter \(r\) is always equal to \(1\). However, when considering isogeometric analysis approximations with polynomial degree \(p\) and regularity \(k\leq p-1\) we have \(r=(p-k)^{d}\)[19; 20]. Notice that a particular case of the previous formula is the case of \(p\) order finite elements in space dimension \(d\) which leads to \(r=p^{d}\)[20; 21], since \(k=0\). Also, the discontinuous Galerkin techniques of degree \(p\) are covered: we have \(r=(p+1)^{d}\)[19] because \(k=-1\).
* The above considerations could be considered also in the case where the fractional Laplacian is defined on a non-Cartesian \(d\)-dimensional domain \(\Omega\), or equipped with variable coefficients, or with approximations on graded grids. In fact, the related GLT theory is already available [3; 18; 19] for encompassing such a generality, while non-Hermitian perturbation tools do not depend on a specific structure of the involved matrix sequences.  
S.S.-C. is responsible for funding acquisition; the contribution of all authors is equal regarding all the listed items: Conceptualization, methodology, validation, investigation, data curation, writing--original draft preparation, writing--review and editing, visualization, supervision, project administration. All authors have read and agreed to the published version of the manuscript.  
The work of [PERSON] is supported by GNCS-INdAM and is funded by the European High-Performance Computing Joint Undertaking (JU) under grant agreement No. 955701. The JU receives support from the European Union's Horizon 2020 research and innovation programme and Belgium, France, Germany, and Switzerland. Furthermore, [PERSON] is grateful for the support of the Laboratory of Theory, Economics and Systems--Department of Computer Science at Athens University of Economics and Business.  
Data Availability Statement: Data are contained within the article.  
\begin{table}
\begin{tabular}{c c c} \hline \hline  & \(\mu(x,y)=1/(x+iy)^{\gamma}\) & \\  & \(\gamma=3\), \(\alpha=1.8\) & \\ \hline \(n\) & - & \(P_{\mathrm{r}}\) \\ \(2^{4}\) & 69 & 17 \\ \(2^{5}\) & 125 & 22 \\ \(2^{6}\) & 236 & 32 \\ \(2^{7}\) & 454 & 47 \\ \(2^{8}\) & \(>\)500 & 70 \\ \(2^{9}\) & \(>\)500 & 108 \\ \hline \hline \end{tabular}
\end{table}
Table 8: Number of preconditioned GMRES iterations to solve the linear system for increasing dimension \(n^{2}\) till \(tol=10^{-11}\).  
**Acknowledgments:** We thank the anonymous Referees for their careful work and for the explicit appreciation of our results.