# Quasinormal Modes in Noncommutative Schwarzschild Black Holes: A Spectral
Analysis

Davide Batic<EMAIL_ADDRESS>Mathematics Department, Khalifa University
of Science and Technology, PO Box 127788, Abu Dhabi, United Arab Emirates
Denys Dutykh<EMAIL_ADDRESS>Mathematics Department, Khalifa University
of Science and Technology, PO Box 127788, Abu Dhabi, United Arab Emirates
Causal Dynamics Pty Ltd, Perth, Australia

###### Abstract

We present a comprehensive analysis of quasinormal modes (QNMs) for
noncommutative geometry-inspired Schwarzschild black holes, encompassing both
non-extreme and extreme cases. By employing a spectral method, we calculate
the QNMs in the context of scalar, electromagnetic, and gravitational
perturbations. Our findings not only challenge previous claims in the
literature regarding the instability of these black holes but also reveal
remarkable stability for both non-extreme and extreme Schwarzschild black
holes under various perturbations.

###### pacs:

04.70.-s,04.70.Bw,04.70.Dy,04.30.-w

##  I Introduction

What exactly is a quasinormal mode (QNM)? In essence, we can think of it as
the real-world counterpart to the idealized normal mode. To illustrate this
concept, let us consider the example of a gently perturbed wine glass. When we
tap the glass lightly, it produces a sound that gradually fades away as energy
dissipates into the surroundings. If this were a normal mode, the ringing
would persist indefinitely. In our current study, we venture into a realm of
objects far more intriguing than an ordinary wine glass, namely, black holes.
When we delve into more complex scenarios, accounting for possible dissipative
phenomena, it becomes natural to expect that an expansion in normal modes is
no longer feasible. In such cases, it will be quite natural for the solutions
to exhibit decay over time instead of oscillating indefinitely. For this
reason, the frequencies will be represented by complex numbers. Let us
consider the possibility that the nature of the equilibrium point of a
physical system around which we perform an expansion is not known a priori.
The QNMs we determine may be characterized by frequencies whose imaginary part
does not describe the damping of solutions around the aforementioned point but
instead signifies an exponential growth over time. It is easy to deduce that,
in this case, the equilibrium point is, in fact, unstable. The significance of
this concept is profound: often, the expansion in QNMs is employed to
ascertain the stability of highly complex systems, and it was precisely this
idea that prompted Regge and Wheeler in their seminal paper Regge and Wheeler
(1957) to initiate the study of black holes within the perturbative framework.

Over the years, various techniques have been developed to study the QNMs, each
with its own peculiarities. One of the most important methods for determining
them is the WKB (Wentzel–Kramers–Brillouin) method Hall (2013). This
particular approach extensively exploits the analogy between the physical
phenomenon we are studying here and the study of energy eigenvalues for wave
functions using the Schrödinger equation in Quantum Mechanics. Another
technique useful for the determination of QNMs makes use of group theory. The
proposed technique, known as the inverted potential method Ferrari and
Mashhoon (1984), is specifically applicable only in the case of a
Schwarzschild black hole, not a Kerr black hole. It is based on the similarity
between the effective and reversed Pöschl–Teller potentials. Furthermore, the
continued fractions method analysis was first introduced by Leaver (1985,
1986) while the Laplace transform method was exploited in Nollert and Schmidt
(1992). In addition to the aforementioned techniques, it is worth mentioning
that the Asymptotic Iterative Method (AIM) has also been applied to determine
the QNMs for Schwarzschild, Reissner-Nordström, and Kerr black holes in the
case of spin $0$, $1/2$, and $2$ perturbations Ciftci et al. (2003); Cho et
al. (2010, 2012); Batic et al. (2018); Macedo (2019); Batic et al. (2019a);
Ciftci et al. (2005). Recently, in a study by Mamani et al. (2022), QNMs of a
Schwarzschild black hole were revisited using an innovative approach called
the spectral method. Spectral methods, though less common in the general
relativity (GR) quasinormal mode (QNM) community, offer substantial advantages
both analytically and numerically. These methods harness spectral expansions
derived from constructive approximation theory. The coefficients within these
expansions are computed numerically, resulting in an approximate analytical
representation of the solution, effectively blending the strengths of both
analytical and numerical approaches.

In contrast to widely used methods such as WKB, inverse potential, and AIM
methods, spectral techniques do not rely on local approximations of the
potential near its maximum nor require constraints on the angular momentum
parameter. The WKB technique also struggles with potentials whose second
derivative at the maximum is zero, thus leading to divergent behaviours in
formulae like (4.1) in Iyer (1987). Finally, for a thorough analysis of the
validity of the WKB approximation in calculating the asymptotic quasinormal
modes of black holes we refer to Daghigh and Green (2012). Asymptotic
iteration methods, while useful, are highly sensitive to the choice of initial
values in the parameter space. An incorrect choice during the initialization
process can affect the stability of the method and the number of iterations
required to compute eigenvalues or quasinormal modes. A recent attempt to find
the optimal selection can be found in Batic and Nowakowski (2023).

On the other hand, spectral methods utilize the full analytical form of the
potential without any need for local approximations or heuristic
modifications. This comprehensive approach not only provides high accuracy but
also ensures exponentially fast convergence for the regular part of the
solution of a radial equation once its singular behaviour has been completely
determined. This feature is particularly advantageous in the search of QNMs of
astrophysical objects. For the aforementioned reasons, in the current paper,
we will employ an improved version of the spectral technique to explore the
QNMs of the Schwarzschild black hole inspired by noncommutative geometry in
the case of a massless particle of spin $s=0$ (scalar perturbation), $s=1$
(electromagnetic perturbation), and $s=2$ (vector-type gravitational
perturbation) Nicolini et al. (2006).

While a substantial body of literature exists dedicated to the investigation
of QNMs in the context of the classic Schwarzschild spacetime Bachelot and
Motet-Bachelot (1993); Gundlach et al. (1994); Andersson (1995, 1997);
Konoplya (2003); Motl and Neitzke (2003); Natário and Schiappa (2004);
Chakrabarti and Gupta (2006); Batic et al. (2018), the situation is
considerably less comprehensive for the noncommutative Schwarzschild manifold.
In the latter case, the literature contains only partial and occasionally
contradictory findings, which we will briefly address in the following
discussion. For example, Giri (2007) analyzed the asymptotic QNMs of the
nonextreme black hole inspired by noncommutative geometry under the assumption
that the black hole mass significantly exceeds its extreme mass. As previously
highlighted by Batic et al. (2019b), in this regime, the QNMs already align
with those of its classical counterpart. Hence, it comes as no surprise that
Giri (2007) concluded that the asymptotic QNMs remain proportional to
$\ln{3}$, thus mirroring the corresponding behaviour in the case of the
classic Schwarzschild solution. Recently, Liang (2018a, b) conducted a study
on QNMs of massless scalar field perturbations in the nonextreme regime
utilizing the third-order WKB approximation. The author specifically focused
on cases with $\ell\in\\{1,2,3\\}$. In this work, the quasinormal frequencies
have been computed for values of the noncommutative parameter $\theta$ across
the interval $[0.01,0.2758]$. If we recall that $\sqrt{\theta}$ represents a
quantum of length, it is reasonable to assume that $\sqrt{\theta}$ is on the
order of the Planck length, $L_{P}$. Consequently, the interval chosen in
Liang (2018a, b) corresponds to a selection of black hole masses within the
range $3.6M_{P}<M<19M_{P}$, where $M_{P}$ represents the Planck mass. In this
specific range and employing the third-order WKB approximation, Liang (2018a,
b) observed stable results for the QNMs. It is noteworthy that, in contrast,
Batic et al. (2019b) identified instabilities using WKB approximations from
third sixth order not only in the narrower range of $1.9M_{P}<M<2.3897M_{P}$
but also in the extreme case $\mu=1.9M_{P}$. Batic et al. (2019b) emphasized
that extending the WKB calculations to the sixth order reveals a lack of
convergence precisely when the presumed instabilities occur. To gain insight
into the nature of these instabilities, Batic et al. (2019b) resorted to
alternative methods. The inverted potential and Asymptotic Iterative Methods
(AIM) failed to detect any instability. However, it is worth noting that both
of these methods entail certain approximations, and the AIM method, in
particular, is known to encounter convergence issues. For a more in-depth
discussion of this aspect, we refer to Batic and Nowakowski (2023). Last but
certainly not least, Batic and Nicolini (2010) conducted an examination of the
stability within the interior of the noncommutative Schwarzschild black hole.
This analysis involved the study of massless scalar field propagation between
the two horizons. Their findings demonstrated that spacetime fuzziness,
influenced by higher momenta of the field, effectively mitigates the classical
exponential blue shift divergence. This mitigation suppresses the emergence of
infinite energy density in regions proximate to the Cauchy horizon.
Conversely, Brown and Mann (2011) investigated the surface gravity of both
inner and outer horizons and posited that the Cauchy horizon might become
unstable when subjected to perturbations from infalling matter. In light of
this ongoing discussion, there is a significant interest in investigating the
(in)stabilities of nonextreme and extreme noncommutative geometry-inspired
black holes by analyzing their respective QNMs. In the present study, we
tackle this challenge by employing the spectral method to analyze quasinormal
frequencies associated with scalar, electromagnetic, and vector-type
gravitational perturbations. Our results unequivocally demonstrate the absence
of instabilities in both the extreme and nonextreme cases. Notably, in the
extreme case, where the event horizon essentially serves as a de facto Cauchy
horizon, our findings align with the conclusions derived in Batic and Nicolini
(2010), further bolstering their validity.

The paper is organised as follows. In Section II, we introduce a suitable
rescaling of the noncommutative Schwarzschild manifold and derive the
effective potential for a massless scalar field in the aforementioned
geometry. Moreover, we also establish the corresponding QNMs boundary
conditions of a massless scalar field in the presence of nonextreme and
extreme noncommutative Schwarzschild black holes. Section III describes the
numerical method used in our study, providing the basis for our computational
approach. In Section IV, we validate our numerical method by applying it to
large-mass noncommutative Schwarzschild black holes. This choice is strategic
as, in this high-mass regime, these black holes converge towards the
characteristics of classical Schwarzschild black holes, providing a robust
benchmark for validation. Here, we also delve into the numerical findings for
both nonextreme and extreme noncommutative black holes. Finally, Section V
offers conclusions and an outlook, reflecting on the implications of our
findings and suggesting directions for future research.

##  II Equations of motion

We consider a massless scalar field $\psi$ immersed in the noncommutative
geometry inspired Schwarzschild background whose line element in units where
$c=G_{N}=1$ reads Nicolini et al. (2006)

$ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\vartheta^{2}+r^{2}\sin^{2}{\vartheta}d\varphi^{2},\quad\vartheta\in[0,\pi],\quad\varphi\in[0,2\pi[,\quad\quad
f(r)=1-\frac{2m(r)}{r}$ (1)

with the mass function

$m(r)=\frac{2M}{\sqrt{\pi}}\gamma\left(\frac{3}{2},\frac{r^{2}}{4\theta}\right),\qquad\gamma\left(\frac{3}{2},\frac{r^{2}}{4\theta}\right)=\int_{0}^{r^{2}/4\theta}dt\sqrt{t}e^{-t},$
(2)

where $M$ is the total mass of a gravitational object, $\theta$ is a parameter
encoding noncommutativity and having the dimension of a length squared. At the
same time, $\gamma(\cdot,\cdot)$ is the incomplete lower gamma function.
Considering the relationship between the upper and lower incomplete gamma
functions as expressed by the formula Abramowitz and Stegun (1972)

$\gamma\left(\frac{3}{2},\frac{r^{2}}{4\theta}\right)+\Gamma\left(\frac{3}{2},\frac{r^{2}}{4\theta}\right)=\frac{\sqrt{\pi}}{2},$
(3)

we can represent the $g_{00}$ metric coefficient in (1) as a combination of
the typical Schwarzschild term augmented by a noncommutativity-induced
perturbation, namely

$f(r)\ =\
1-\frac{2M}{r}+\frac{4M}{\sqrt{\pi}r}\>\Gamma\left(\frac{3}{2},\frac{r^{2}}{4\theta}\right).$
(4)

Note that when $r/\sqrt{\theta}\to+\infty$, the term
$\Gamma(3/2,r^{2}/4\theta)$ tends to zero, yielding the classic Schwarzschild
metric. If we introduce the following rescaling

$x=\frac{r}{2M},\qquad\mu=\frac{M}{\sqrt{\theta}}$ (5)

coupled with the identities Abramowitz and Stegun (1972)

$\gamma\left(\frac{3}{2},w^{2}\right)=\frac{1}{2}\gamma\left(\frac{1}{2},w^{2}\right)-we^{-w^{2}},\quad\gamma\left(\frac{1}{2},w^{2}\right)=\sqrt{\pi}\mbox{erf}(w),$
(6)

we can cast (4) in the form

$f(x)\ =\ 1-\frac{\mbox{erf}\left(\mu
x\right)}{x}+\frac{2\mu}{\sqrt{\pi}}e^{-\mu^{2}x^{2}},$ (7)

where $\mbox{erf}(\cdot)$ represents the error function. The picture emerging
from Figure 1 where $f$ has been plotted against $x$, is the following: an
extreme noncommutative geometry inspired Schwarzschild black hole appears when
the rescaled mass takes the critical value $\mu_{e}=1.904119076\ldots$. In
this regime, the Cauchy horizon $x_{c}$ and the event horizon $x_{h}$
coincide, and their values correspond to
$x_{e}=x_{c}=x_{h}=0.7936575898\ldots$ Whenever $\mu<\mu_{e}$, the line
element (1) represents a naked gravitational droplet regular at $x=0$ Nicolini
et al. (2006). Finally, two distinct horizons exist for $\mu>\mu_{e}$.
Differently as in the case of the Schwarzschild metric, if we probe into
shorter distances, _i.e._ $r\ll\sqrt{\theta}$, it is possible to show the
absence of a curvature central singularity which is instead replaced by a
regular deSitter core Nicolini et al. (2006). In the present work, we focus
our analysis on the QNMs for the nonextreme and extreme cases.

Figure 1:  Plot of the metric coefficient $g_{00}=f(x)$ defined in (7). An
extreme black hole occurs when $\mu_{e}=1.904119076\ldots$ (solid line) and
the horizon is located at $x_{e}=0.7936575898\ldots$ A nonextreme black hole
with two horizons is present for $\mu>\mu_{e}$ (see the dotted line
corresponding to the choice $\mu=2.3$). A naked gravitational droplet occurs
when $0<\mu<\mu_{e}$ (see the dashed line for $\mu=1.3$).

Considering this, it is advantageous to implement the rescaling $z=x/x_{h}$,
as this approach effectively maps the event horizon to $1$ in both the extreme
and nonextreme cases. Hence, the $g_{00}$ metric coefficient, as expressed by
(7), becomes

$f(z)=1-\frac{\mbox{erf}\left(\mu
x_{h}z\right)}{x_{h}z}+\frac{2\mu}{\sqrt{\pi}}e^{-\mu^{2}x_{h}^{2}z^{2}},$ (8)

while in the presence of a noncommutative geometry-inspired Schwarzschild
black hole, the equation governing a massless Klein–Gordon field (assumed to
have a time dependence of the form $e^{-i\omega t}$ and an angular component
described by spherical harmonics) is as follows Batic et al. (2019b)

$f(r)\frac{d}{dr}\left(f(r)\frac{d\psi_{\omega\ell\epsilon}}{dr}\right)+\left[\omega^{2}-U_{\epsilon}(r)\right]\psi_{\omega\ell\epsilon}(r)=0,\qquad
U_{\epsilon}(r)=f(r)\left[\frac{\epsilon}{r}\frac{df}{dr}+\frac{\ell(\ell+1)}{r^{2}}\right],\qquad\epsilon=1-s^{2}$
(9)

with $\ell=0,1,2,\ldots$ and $\epsilon=-1$ (massless scalar perturbation
$s=0$), $\epsilon=0$ (electromagnetic perturbation $s=1$), and $\epsilon=3$
(vector-type gravitational perturbation $s=2$). A comprehensive analysis of
the effective potential, along with its graphical representations for both
configurations, has been extensively covered in Batic et al. (2019b). By means
of the substitution $r=2Mx_{h}z$, the above equation can be recast in the
equivalent form

$\frac{f(z)}{4}\frac{d}{dz}\left(f(z)\frac{d\psi_{\Omega\ell\epsilon}}{dz}\right)+\left[x_{h}^{2}\Omega^{2}-V_{\epsilon}(z)\right]\psi_{\Omega\ell\epsilon}(z)=0,\qquad
V_{\epsilon}(z)=\frac{f(z)}{4}\left[\frac{\epsilon}{z}\frac{df}{dz}+\frac{\ell(\ell+1)}{z^{2}}\right],\qquad\Omega=M\omega$
(10)

with $f(z)$ as given in (8). As a consistency check, note that in the limit
corresponding to the classic Schwarzschild case, _i.e._ $x_{h}\to 1$ and
$r/\sqrt{\theta}\to\infty$, the above equation correctly reproduces equation
$(2)$ in Leaver (1985), where the term $\omega^{3}$ should be corrected to
$\omega^{2}$. In the following analysis, we focus on computing the QNMs for
the spectral problem stated in (10). For this purpose, we represent $\Omega$
as $\Omega=\Omega_{R}+i\Omega_{I}$, where $\Omega_{I}<0$ ensures that the
perturbation is dumped in time. The boundary conditions are set so that the
radial field exhibits inward radiation at the event horizon and outward
radiation at spatial infinity. This necessitates a detailed examination of the
solution asymptotic behaviour in (10), both near the event horizon ($z\to
1^{+}$) and at large spatial distances ($z\to+\infty$). Moreover, to compute
the QNMs using the spectral method, we must recast the differential equation
in (10) and the appropriate boundary conditions over the compact interval
$[-1,1]$. This adjustment is essential as the method expands the regular part
of the eigenfunctions by means of Chebyshev polynomials.

###  II.1 The non-extreme case

This scenario focuses on mass parameters $\mu>\mu_{e}$. The metric coefficient
$g_{00}$ is given as defined in (8). It is important to note that the Cauchy
and event horizons are distinct in the present case, and therefore, $g_{00}$
exhibits a simple zero at $z=1$. To establish the QNM boundary conditions at
the event horizon and at infinity, we first need to determine the asymptotic
behaviour of the radial solution $\psi_{\Omega\ell\epsilon}$ as $z\to 1^{+}$
and as $z\to+\infty$. We can then extract the QNM boundary conditions from
this asymptotic data. As a consistency check, we will also verify that in the
limit $r/\sqrt{\theta}$, these conditions correctly reproduce the
corresponding ones for the classic Schwarzschild case specified in equation
$(4)$ of Leaver (1985). We split our analysis by examining the behaviour of
the radial field in two different regions.

  1. 1.

Asymptotic behaviour as $z\to 1^{+}$: Given that $z=1$ is a simple zero of
$f(z)$, we can represent the latter in the form $f(z)=(z-1)g(z)$ where $g$ is
an analytic function at $z=1$ with the property that $g(1)=f^{\prime}(1)\neq
0$. Here, the prime symbol stands for differentiation with respect to $z$.
This representation enables us to reformulate (10) in the form

$\displaystyle\frac{d^{2}\psi_{\Omega\ell\epsilon}}{dz^{2}}+p(z)\frac{d\psi_{\Omega\ell\epsilon}}{dz}+q(z)\psi_{\Omega\ell\epsilon}(z)=0,$
(11) $\displaystyle p(z)=\frac{f^{{}^{\prime}}(z)}{(z-1)g(z)},\quad
q(z)=\frac{4x_{h}^{2}\Omega^{2}}{(z-1)^{2}g^{2}(z)}-\frac{1}{(z-1)g(z)}\left[\frac{\epsilon}{z}f^{{}^{\prime}}(z)+\frac{\ell(\ell+1)}{z^{2}}\right].$
(12)

Since $p$ and $q$ have poles of order one and two at $z=1$, respectively, this
point is classified as a regular singular point of (11), according to
Frobenius theory Ince (1956). Hence, we can construct solutions of the form

$\psi_{\Omega\ell\epsilon}(z)=(z-1)^{\rho}\sum_{\kappa=0}^{\infty}a_{\kappa}(z-1)^{\kappa}.$
(13)

The leading behavior at $z=1$ is represented by the term $(z-1)^{\rho}$ where
$\rho$ is determined by the indicial equation

$\rho(\rho-1)+P_{0}\rho+Q_{0}=0$ (14)

with

$P_{0}=\lim_{z\to 1}(z-1)p(z)=1,\qquad Q_{0}=\lim_{z\to
1}(z-1)^{2}q(z)=\left(\frac{2x_{h}\Omega}{f^{{}^{\prime}}(1)}\right)^{2}.$
(15)

The roots of (14) are $\rho_{\pm}=\pm 2ix_{h}\Omega/f^{{}^{\prime}}(1)$ and
the correct QNM boundary condition at $z=1$ reads

$\psi_{\Omega\ell\epsilon}\underset{{z\to
1^{+}}}{\longrightarrow}(z-1)^{-2ix_{h}\alpha\Omega},\quad\alpha=\frac{1}{f^{{}^{\prime}}(1)}.$
(16)

We can easily obtain a compact formula for $f^{{}^{\prime}}(1)$ if we observe
that

$f^{{}^{\prime}}(1)=\frac{\mbox{erf}(\mu
x_{h})}{x_{h}}-\frac{2\mu}{\sqrt{\pi}}(1+2\mu^{2}x_{h}^{2})e^{-\mu^{2}x_{h}^{2}}$
(17)

can be further simplified by taking into account the condition $f(1)=0$.
Solving the latter with respect to the error function leads to

$\mbox{erf}(\mu
x_{h})=x_{h}\left(1+\frac{2\mu}{\sqrt{\pi}}e^{-\mu^{2}x_{h}^{2}}\right),$ (18)

which in turns allows to express (17) as follows

$f^{{}^{\prime}}(1)=1-\frac{4\mu^{3}x_{h}^{2}}{\sqrt{\pi}}e^{-\mu^{2}x_{h}^{2}}.$
(19)

At this point, a remark is in order. We note that as $\sqrt{\theta}\to 0$, the
rescaled mass parameter $\mu=M/\sqrt{\theta}\to+\infty$. This implies that
$f^{\prime}(1)\to 1$ aligns with the expected behaviour in the classic
Schwarzschild case. Applying this observation to (16) and considering that in
this limit $x_{h}\to 1$, it becomes evident that the prescription (16)
accurately reflects the QNM boundary condition at the event horizon, as
delineated in equation $(4)$ of Leaver (1985)

  2. 2.

Asymptotic behaviour as $z\to+\infty$: We first rewrite (10) as

$\frac{d^{2}\psi_{\Omega\ell\epsilon}}{dz^{2}}+P(z)\frac{d\psi_{\Omega\ell\epsilon}}{dz}+Q(z)\psi_{\Omega\ell\epsilon}(z)=0,\qquad
P(z)=\frac{f^{{}^{\prime}}(z)}{f(z)},\qquad
Q(z)=\frac{4\left[x_{h}^{2}\Omega^{2}-V_{\epsilon}(z)\right]}{f^{2}(z)}$ (20)

with $V_{\epsilon}(z)$ defined in (10). In order to classify the point at
infinity, it is convenient to introduce the transformation $\tau=1/z$. Then,
the above differential equation becomes

$\frac{d^{2}\psi_{\Omega\ell\epsilon}}{d\tau^{2}}+C(\tau)\frac{d\psi_{\Omega\ell\epsilon}}{d\tau}+D(\tau)\psi_{\Omega\ell\epsilon}(\tau)=0,\quad
C(\tau)=\frac{2}{\tau}+\frac{1}{f(\tau)}\frac{df}{d\tau},\quad
D(\tau)=\frac{4\left[x_{h}^{2}\Omega^{2}-V_{\epsilon}(\tau)\right]}{\tau^{4}f^{2}(\tau)}.$
(21)

Since for $k=1$ we have

$\displaystyle\tau^{k+1}C(\tau)$ $\displaystyle=$ $\displaystyle
C_{0}+C_{1}\tau+\mathcal{O}(\tau^{2}),\qquad C_{0}=0,\qquad C_{1}=2,$ (22)
$\displaystyle\tau^{2k+2}D(\tau)$ $\displaystyle=$ $\displaystyle
D_{0}+D_{1}\tau+\mathcal{O}(\tau^{2}),\qquad D_{0}=4x_{h}^{2}\Omega^{2},\quad
D_{1}=8x_{h}\Omega^{2}$ (23)

with $D_{0}\neq 0$, according to Bender and Orszag (1999), the point at
infinity is an irregular singular point of rank one. Consequently, the
asymptotic behaviour of the solutions to equation (20) can be efficiently
deduced using the method outlined in Olver (1994). For this purpose, we start
by observing that

$P(z)=\sum_{\kappa=0}^{\infty}\frac{\mathfrak{f}_{\kappa}}{z^{k}}=\mathcal{O}\left(\frac{1}{z^{2}}\right),\qquad
Q(z)=\sum_{\kappa=0}^{\infty}\frac{\mathfrak{g}_{\kappa}}{z^{k}}=4x_{h}^{2}\Omega^{2}+\frac{8x_{h}\Omega^{2}}{z}+\mathcal{O}\left(\frac{1}{z^{2}}\right).$
(24)

Given that at least one of the coefficients $\mathfrak{f}_{0}$,
$\mathfrak{g}_{0}$, $\mathfrak{g}_{1}$ is nonzero, a formal solution to (20)
is represented by Olver (1994)

$\psi^{(j)}_{\Omega\ell\epsilon}(z)=z^{\mu_{j}}e^{\lambda_{j}z}\sum_{\kappa=0}^{\infty}\frac{a_{\kappa,j}}{z^{\kappa}},\qquad
j\in\\{1,2\\},$ (25)

where $\lambda_{1}$, $\lambda_{2}$, $\mu_{1}$ and $\mu_{2}$ are the roots of
the characteristic equations

$\lambda^{2}+\mathfrak{f}_{0}\lambda+\mathfrak{g}_{0}=0,\quad\mu_{j}=-\frac{\mathfrak{f}_{1}\lambda_{j}+\mathfrak{g}_{1}}{\mathfrak{f}_{0}+2\lambda_{j}}.$
(26)

A straightforward computation shows that $\lambda_{\pm}=\pm 2ix_{h}\Omega$ and
$\mu_{\pm}=\pm 2i\Omega$. As a result, the QNM boundary condition at space-
like infinity can be expressed, as

$\psi_{\Omega\ell\epsilon}\underset{{z\to+\infty}}{\longrightarrow}z^{2i\Omega}e^{2ix_{h}\Omega
z}.$ (27)

It is gratifying to observe that it reproduces correctly the corresponding
condition for the classic Schwarzschild case.

At this point, we can proceed as in Leaver (1985), that is, we transform the
radial function $\psi_{\Omega\ell\epsilon}(z)$ into a new radial function
$\Phi_{\Omega\ell\epsilon}(z)$ such that the QNM boundary conditions are
automatically implemented and $\Phi_{\sigma\ell\epsilon}(z)$ is regular at
$z=1$ and at space-like infinity. To this aim, we consider the transformation

$\psi_{\Omega\ell\epsilon}(z)=z^{2i(1+x_{h}\alpha)\Omega}(z-1)^{-2ix_{h}\alpha\Omega}e^{2ix_{h}\Omega(z-1)}\Phi_{\Omega\ell\epsilon}(z),\quad\alpha=\frac{1}{1-\frac{4\mu^{3}x_{h}^{2}}{\sqrt{\pi}}e^{-\mu^{2}x_{h}^{2}}}.$
(28)

Note that in the limit of a vanishing noncommutative parameter, the above
ansatz reduces to the corresponding one represented by equation (5) in Leaver
(1985) for the classic Schwarzschild case. If we substitute (28) into (10), we
end up with the following ordinary differential equation for the radial
eigenfunctions, namely

$P_{2}(z)\Phi^{{}^{\prime\prime}}_{\Omega\ell\epsilon}(z)+P_{1}(z)\Phi^{{}^{\prime}}_{\Omega\ell\epsilon}(z)+P_{0}(z)\Phi_{\Omega\ell\epsilon}(z)=0$
(29)

with

$\displaystyle P_{2}(z)$ $\displaystyle=$
$\displaystyle\frac{z^{2}(z-1)^{2}}{4}f^{2}(z),$ (30) $\displaystyle P_{1}(z)$
$\displaystyle=$ $\displaystyle
z(z-1)f(z)\left\\{\frac{z(z-1)}{4}f^{{}^{\prime}}(z)-i\Omega
f(z)\left[1+x_{h}\alpha+(x_{h}-1)z-x_{h}z^{2}\right]\right\\},$ (31)
$\displaystyle P_{0}(z)$ $\displaystyle=$
$\displaystyle-\Omega^{2}Q_{+}(z)Q_{-}(z)+i\Omega
f(z)L(z)-z^{2}(z-1)^{2}V_{\epsilon}(z),$ (32) $\displaystyle Q_{\pm}(z)$
$\displaystyle=$ $\displaystyle x_{h}z(z-1)\left[f(z)\pm
1\right]-f(z)\left(1+x_{h}\alpha-z\right),$ (33) $\displaystyle L(z)$
$\displaystyle=$
$\displaystyle\frac{(z-1)^{2}}{2}\left[z(1+x_{h}z)f^{{}^{\prime}}(z)-f(z)\right]-\frac{x_{h}\alpha}{2}\left[z(z-1)f^{{}^{\prime}}(z)+(1-2z)f(z)\right].$
(34)

Let us now introduce the transformation $z=2/(1-y)$ mapping the point at
infinity and the event horizon to $y=1$ and $y=-1$, respectively. Furthermore,
a dot denotes differentiation with respect to the new variable $y$. Then,
equation (29) becomes

$S_{2}(y)\ddot{\Phi}_{\Omega\ell\epsilon}(y)+S_{1}(y)\dot{\Phi}_{\Omega\ell\epsilon}(y)+S_{0}(y)\Phi_{\Omega\ell\epsilon}(y)=0,$
(35)

where

$\displaystyle S_{2}(y)$ $\displaystyle=$
$\displaystyle\frac{(1+y)^{2}}{4}f^{2}(y),$ (36) $\displaystyle S_{1}(y)$
$\displaystyle=$ $\displaystyle
i\Omega\frac{1+y}{(1-y)^{2}}f^{2}(y)\left[(1+y)(1+2x_{h}-y)-x_{h}\alpha(1-y)^{2}\right]-\frac{(1+y)^{2}}{2(1-y)}f^{2}(y)+\frac{(1+y)^{2}}{4}f(y)\dot{f}(y),$
(37) $\displaystyle S_{0}(y)$ $\displaystyle=$
$\displaystyle\Omega^{2}\Sigma_{2}(y)+i\Omega\Sigma_{1}(y)+\Sigma_{0}(y)$ (38)

with

$\displaystyle\Sigma_{2}(y)$ $\displaystyle=$
$\displaystyle\frac{4x_{h}^{2}(1+y)^{2}}{(1-y)^{4}}-\frac{f^{2}(y)}{(1-y)^{4}}\left[(1+y)(1+2x_{h}-y)-x_{h}\alpha(1-y)^{2}\right]^{2},$
(39) $\displaystyle\Sigma_{1}(y)$ $\displaystyle=$
$\displaystyle\frac{f(y)}{2}\left\\{\left(\frac{1+y}{1-y}\right)^{2}\left[(1+2x_{h}-y)\dot{f}(y)-f(y)\right]+\frac{x_{h}\alpha}{1-y}\left[(3+y)f(y)-(1-y^{2})\dot{f}(y)\right]\right\\},$
(40) $\displaystyle\Sigma_{0}(y)$ $\displaystyle=$
$\displaystyle-\frac{4(1+y)^{2}}{(1-y)^{4}}V_{\epsilon}(y).$ (41)

Notice that we must also require that $\Phi_{\Omega\ell\epsilon}(y)$ is
regular at $y=\pm 1$. As a result of the transformation introduced above, we
have

$f(y)=1-\frac{1-y}{2x_{h}}\mbox{erf}\left(\frac{2\mu
x_{h}}{1-y}\right)+\frac{2\mu}{\sqrt{\pi}}e^{-\frac{4\mu^{2}x_{h}^{2}}{(1-y)^{2}}},\qquad
V_{\epsilon}(y)=\frac{(1-y)^{2}}{16}f(y)\left[\epsilon(1-y)\dot{f}(y)+\ell(\ell+1)\right].$
(42)

Table 1: Classification of the points $y=\pm 1$ for the relevant functions defined by (42), (36), (37) and (38). The abbreviations $z$ ord $n$ and $p$ ord $m$ stand for zero of order $n$ and pole of order $m$, respectively. $y$ | $f(y)$ | $V_{\epsilon}(y)$ | $S_{2}(y)$ | $S_{1}(y)$ | $S_{0}(y)$
---|---|---|---|---|---
$-1$ | z ord 1 | z ord 1 | z ord 4 | z ord 3 | z ord 3
$+1$ | $+1$ | z ord 2 | $+1$ | p ord 2 | p ord 2

Table 1 shows that the coefficients of the differential equation (35) share a
common zero of order $3$ at $y=-1$ while $y=1$ is a pole of order $2$ for the
coefficients $S_{1}$ and $S_{0}$. Hence, in order to apply the spectral
method, we need to multiply (35) by $(1-y)^{2}/(1+y)^{3}$. As a result, we end
up with the following differential equation

$M_{2}(y)\ddot{\Phi}_{\Omega\ell\epsilon}(y)+M_{1}(y)\dot{\Phi}_{\Omega\ell\epsilon}(y)+M_{0}(y)\Phi_{\Omega\ell\epsilon}(y)=0,$
(43)

where

$M_{2}(y)=\frac{(1-y)^{2}}{4(1+y)}f^{2}(y),\qquad M_{1}(y)=i\Omega
N_{1}(y)+N_{0}(y),\qquad M_{0}(y)=\Omega^{2}C_{2}(y)+i\Omega
C_{1}(y)+C_{0}(y)$ (44)

with

$\displaystyle N_{1}(y)$ $\displaystyle=$ $\displaystyle
f^{2}(y)\left[\frac{1+2x_{h}-y}{1+y}-x_{h}\alpha\left(\frac{1-y}{1+y}\right)^{2}\right],\quad
N_{0}(y)=\frac{f(y)}{4(1+y)}\frac{d}{dy}\left((1-y)^{2}f(y)\right),$ (45)
$\displaystyle C_{2}(y)$ $\displaystyle=$
$\displaystyle\frac{4x_{h}^{2}}{(1+y)(1-y)^{2}}-\frac{f^{2}(y)}{(1+y)^{3}(1-y)^{2}}\left[(1+y)(1+2x_{h}-y)-x_{h}\alpha(1-y)^{2}\right]^{2},$
(46) $\displaystyle C_{1}(y)$ $\displaystyle=$
$\displaystyle\frac{f(y)}{2(1+y)}\left\\{(1+2x_{h}-y)\dot{f}(y)-f(y)+x_{h}\alpha\frac{1-y}{(1+y)^{2}}\left[(3+y)f(y)-(1-y^{2})\dot{f}(y)\right]\right\\},$
(47) $\displaystyle C_{0}(y)$ $\displaystyle=$
$\displaystyle-\frac{4V_{\epsilon}(y)}{(1+y)(1-y)^{2}}.$ (48)

It can be easily verified with Maple that

$\displaystyle\lim_{y\to 1^{-}}M_{2}(y)=0=\lim_{y\to-1^{+}}M_{2}(y),$ (49)
$\displaystyle\lim_{y\to
1^{-}}M_{1}(y)=ix_{h}\Omega,\quad\lim_{y\to-1^{+}}M_{1}(y)=i\Omega\Lambda_{1}+\Lambda_{0},$
(50) $\displaystyle\lim_{y\to
1^{-}}M_{0}(y)=A_{2}\Omega^{2}+A_{0},\quad\lim_{y\to-1^{+}}M_{0}(y)=B_{2}\Omega^{2}+i\Omega
B_{1}+B_{0},$ (51)

where

$\displaystyle\Lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{x_{h}\left(4\mu^{3}x_{h}^{2}e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi}\right)}{\sqrt{\pi}},\quad\Lambda_{0}=\frac{\left(4\mu^{3}x_{h}^{2}e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi}\right)^{2}}{4\pi},$
(52) $\displaystyle A_{2}$ $\displaystyle=$ $\displaystyle
1-\frac{\sqrt{\pi}x_{h}^{2}}{4\mu^{3}x_{h}^{2}e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi}},\quad
A_{0}=-\frac{\ell(\ell+1)}{8},$ (53) $\displaystyle B_{2}$ $\displaystyle=$
$\displaystyle-\frac{x_{h}(1+x_{h})}{\sqrt{\pi}}(4\mu^{3}x_{h}^{2}e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi})+\frac{2x_{h}^{2}\left[2\mu^{3}x_{h}^{2}(1+\mu^{2}x_{h}^{2})e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi}\right]}{4\mu^{3}x_{h}^{2}e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi}},$
(54) $\displaystyle B_{1}$ $\displaystyle=$
$\displaystyle\frac{\mathfrak{a}_{3}e^{-3\mu^{2}x_{h}^{2}}+\mathfrak{a}_{2}e^{-2\mu^{2}x_{h}^{2}}+\mathfrak{a}_{1}e^{-\mu^{2}x_{h}^{2}}+\mathfrak{a}_{0}}{4\pi(4\mu^{3}x_{h}^{2}e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi})},$
(55) $\displaystyle B_{0}$ $\displaystyle=$
$\displaystyle\frac{4\mu^{3}x_{h}^{2}e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi}}{8\sqrt{\pi}}\left[\ell(\ell+1)-\frac{\epsilon}{\sqrt{\pi}}(4\mu^{3}x_{h}^{2}e^{-\mu^{2}x_{h}^{2}}-\sqrt{\pi})\right]$
(56)

and

$\displaystyle\mathfrak{a}_{3}$ $\displaystyle=$ $\displaystyle
64\mu^{9}x_{h}^{6}(1+x_{h}),\quad\mathfrak{a}_{2}=-16\sqrt{\pi}\mu^{6}x_{h}^{4}(\mu^{2}x_{h}^{3}+4x_{h}+3),$
(57) $\displaystyle\mathfrak{a}_{1}$ $\displaystyle=$ $\displaystyle
4\pi\mu^{3}x_{h}^{2}(\mu^{2}x_{h}^{3}+6x_{h}+3),\quad\mathfrak{a}_{0}=-\pi^{3/2}(1+3x_{h}).$
(58)

In the final step leading to the application of the spectral method, we recast
the differential equation (43) into the following form

$\widehat{L}_{0}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]+i\widehat{L}_{1}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]\Omega+\widehat{L}_{2}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]\Omega^{2}=0$
(59)

with

$\displaystyle\widehat{L}_{0}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]$
$\displaystyle=$
$\displaystyle\widehat{L}_{00}(y)\Phi_{\Omega\ell\epsilon}+\widehat{L}_{01}(y)\dot{\Phi}_{\Omega\ell\epsilon}+\widehat{L}_{02}(y)\ddot{\Phi}_{\Omega\ell\epsilon},$
(60)
$\displaystyle\widehat{L}_{1}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]$
$\displaystyle=$
$\displaystyle\widehat{L}_{10}(y)\Phi_{\Omega\ell\epsilon}+\widehat{L}_{11}(y)\dot{\Phi}_{\Omega\ell\epsilon}+\widehat{L}_{12}(y)\ddot{\Phi}_{\Omega\ell\epsilon},$
(61)
$\displaystyle\widehat{L}_{2}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]$
$\displaystyle=$
$\displaystyle\widehat{L}_{20}(y)\Phi_{\Omega\ell\epsilon}+\widehat{L}_{21}(y)\dot{\Phi}_{\Omega\ell\epsilon}+\widehat{L}_{22}(y)\ddot{\Phi}_{\Omega\ell\epsilon}.$
(62)

Moreover, in Table 2, we have summarized the $\widehat{L}_{ij}$ appearing in
(60)–(62) and their limiting values at $y=\pm 1$.

Table 2: Definitions of the coefficients $\widehat{L}_{ij}$ and their corresponding behaviours at the endpoints of the interval $-1\leq y\leq 1$. The symbols appearing in this table have been defined in (44)-(56). $(i,j)$ | $\displaystyle{\lim_{y\to-1^{+}}}\widehat{L}_{ij}$ | $\widehat{L}_{ij}$ | $\displaystyle{\lim_{y\to 1^{-}}}\widehat{L}_{ij}$
---|---|---|---
$(0,0)$ | $B_{0}$ | $C_{0}$ | $A_{0}$
$(0,1)$ | $\Lambda_{0}$ | $N_{0}$ | $0$
$(0,2)$ | $0$ | $M_{2}$ | $0$
$(1,0)$ | $B_{1}$ | $C_{1}$ | $0$
$(1,1)$ | $\Lambda_{1}$ | $N_{1}$ | $x_{h}$
$(1,2)$ | $0$ | $0$ | $0$
$(2,0)$ | $B_{2}$ | $C_{2}$ | $A_{2}$
$(2,1)$ | $0$ | $0$ | $0$
$(2,2)$ | $0$ | $0$ | $0$

###  II.2 The extreme case

In this scenario characterized by $\mu=\mu_{e}$, we still introduce the
rescaling $x=r/(2M)$. However, since the Cauchy and the event horizon
coincide, _i.e._ $x_{e}=x_{c}=x_{h}=0.7936575898\ldots$, it is useful to
introduce a second rescaling given by $\xi=x/x_{e}$. Hence, the $g_{00}$
metric coefficient, as expressed by (7), becomes

$f_{e}(\xi)=1-\frac{\mbox{erf}\left(\mu_{e}x_{e}\xi\right)}{x_{e}\xi}+\frac{2\mu_{e}}{\sqrt{\pi}}e^{-\mu_{e}^{2}x_{e}^{2}\xi^{2}},$
(63)

while the radial part of a massless Klein–Gordon field can be expressed in the
equivalent form

$\frac{f_{e}(\xi)}{4}\frac{d}{d\xi}\left(f_{e}(\xi)\frac{d\psi_{\Omega\ell\epsilon}}{d\xi}\right)+\left[x_{e}^{2}\Omega^{2}-V_{\epsilon}(\xi)\right]\psi_{\Omega\ell\epsilon}(\xi)=0,\quad
V_{\epsilon}(\xi)=\frac{f_{e}(\xi)}{4}\left[\frac{\epsilon}{\xi}\frac{df_{e}}{d\xi}+\frac{\ell(\ell+1)}{\xi^{2}}\right],\quad\Omega=M\omega$
(64)

with $f_{e}(\xi)$ as given in (63). The fact that $f_{e}(\xi)$ exhibits a zero
of order two at $\xi=1$, implies that $f_{e}(1)=0=f^{{}^{\prime}}(1)$ where
the prime denotes differentiation with respect to $\xi$. This observation
allows us to derive the following functional relations

$e^{-\mu_{e}^{2}x_{e}^{2}}=\frac{\sqrt{\pi}}{4\mu_{e}^{3}x_{e}^{2}},\quad\mbox{erf}\left(\mu_{e}x_{e}\right)=\frac{1+2\mu_{e}^{2}x_{e}^{2}}{2\mu_{e}^{2}x_{e}},$
(65)

which, in turn, plays an important role in simplifying the forthcoming
computations. In order to derive the QNM boundary conditions at the event
horizon and at infinity, we first need to determine the asymptotic behaviour
of the radial solution $\psi_{\Omega\ell\epsilon}$ as $\xi\to 1^{+}$ and as
$\xi\to+\infty$. We can then accurately extract the QNM boundary conditions
from this asymptotic data.

  1. 1.

Asymptotic behaviour as $\xi\to 1^{+}$: Taking into account that $\xi=1$ is a
double zero of $f_{e}(\xi)$, we can represent the latter in the form
$f_{e}(\xi)=(\xi-1)^{2}h(\xi)$ where $h$ is an analytic function at $\xi=1$
with the property that
$h(1)=f_{e}^{{}^{\prime\prime}}(1)/2=\mu_{e}^{2}x_{e}^{2}-1\approx 1.2837$.
This representation enables us to reformulate (64) in the form

$\displaystyle\frac{d^{2}\psi_{\Omega\ell\epsilon}}{d\xi^{2}}+\mathfrak{p}(\xi)\frac{d\psi_{\Omega\ell\epsilon}}{d\xi}+\mathfrak{q}(\xi)\psi_{\Omega\ell\epsilon}(\xi)=0,$
(66)
$\displaystyle\mathfrak{p}(\xi)=\frac{2}{\xi-1}+\frac{h^{{}^{\prime}}(\xi)}{h(\xi)},\quad\mathfrak{q}(\xi)=\frac{4x_{e}^{2}\Omega^{2}}{(\xi-1)^{4}h^{2}(\xi)}-\frac{1}{(\xi-1)^{2}h(\xi)}\left[\frac{\epsilon}{\xi}f_{e}^{{}^{\prime}}(\xi)+\frac{\ell(\ell+1)}{\xi^{2}}\right].$
(67)

Because $\mathfrak{q}$ has a fourth-order pole at $\xi=1$, it follows that
$\xi=1$ qualifies as an irregular singularity, rendering Frobenius’s theory
inapplicable in this scenario. On the other hand, since for $k=1$ we have

$(\xi-1)^{k+1}\mathfrak{p}(\xi)=\mathcal{O}(\xi-1),\quad(\xi-1)^{2k+2}\mathfrak{q}(\xi)=\mathfrak{d}_{0}+\mathcal{O}(\xi-1)^{2},\quad\mathfrak{d}_{0}=\frac{4x_{e}^{2}\Omega^{2}}{h^{2}(1)}$
(68)

with $\mathfrak{d}_{0}\neq 0$, then, according to Bender and Orszag (1999),
$\xi=1$ is an irregular singular point of rank one. Consequently, the leading
behaviour of the solutions to equation (66) in a neighbourhood of the event
horizon can be efficiently deduced using the method outlined in Olver (1994).
To this purpose, we start by observing that by means of the transformation
$\tau=(\xi-1)^{-1}$ mapping the event horizon at infinity and infinity to
zero, (66) becomes

$\displaystyle\frac{d^{2}\psi_{\Omega\ell\epsilon}}{d\tau^{2}}+\mathfrak{C}(\tau)\frac{d\psi_{\Omega\ell\epsilon}}{d\tau}+\mathfrak{D}(\tau)\psi_{\Omega\ell\epsilon}(\tau)=0,$
(69)
$\displaystyle\mathfrak{C}(\tau)=\sum_{\kappa=0}^{\infty}\frac{\mathfrak{c}_{\kappa}}{\tau^{k}}=\mathcal{O}\left(\frac{1}{\tau^{2}}\right),\quad\mathfrak{D}(\tau)=\sum_{\kappa=0}^{\infty}\frac{\mathfrak{d}_{\kappa}}{\tau^{k}}=\frac{4x_{e}^{2}\Omega^{2}}{(\mu_{e}^{2}x_{e}^{2}-1)^{2}}+\frac{16x_{e}^{2}\Omega^{2}(\mu_{e}^{4}x_{e}^{4}-\mu_{e}^{2}x_{e}^{2}-1)}{3(\mu_{e}^{2}x_{e}^{2}-1)^{3}\tau}+\mathcal{O}\left(\frac{1}{\tau^{2}}\right).$
(70)

Since at least one of the coefficients $\mathfrak{c}_{0}$, $\mathfrak{d}_{0}$,
$\mathfrak{d}_{1}$ is nonzero, a formal solution to (66) is given by Olver
(1994)

$\psi^{(\pm)}_{\Omega\ell\epsilon}(\tau)=\tau^{\mu_{\pm}}e^{\lambda_{\pm}\tau}\sum_{\kappa=0}^{\infty}\frac{\mathfrak{a}_{\kappa,\pm}}{\tau^{\kappa}},$
(71)

where $\lambda_{\pm}$, and $\mu_{\pm}$ are the roots of the characteristic
equations

$\lambda_{\pm}^{2}+\mathfrak{d}_{0}=0,\quad\mu_{\pm}=-\frac{\mathfrak{d}_{1}}{2\lambda_{\pm}}.$
(72)

A straightforward computation shows that

$\lambda_{\pm}=\pm\frac{2ix_{e}\Omega}{\mu_{e}^{2}x_{e}^{2}-1},\quad\mu_{\pm}=\pm\frac{4ix_{e}\Omega(\mu_{e}^{4}x_{e}^{4}-\mu_{e}^{2}x_{e}^{2}-1)}{3(\mu_{e}^{2}x_{e}^{2}-1)^{2}}.$
(73)

At this point, it is important to note that a radial field exhibiting solely
inward radiation near the event horizon ($\xi\to 1^{+}$) corresponds, under
the transformation $\tau=(\xi-1)^{-1}$, to an outward radiating field as we
approach spatial infinity ($\tau\to+\infty^{-}$). This indicates that we
should choose the plus sign in the formulas above. Consequently, the correct
QNM boundary condition at $\xi=1$ is as follows

$\psi_{\Omega\ell\epsilon}\underset{{\xi\to
1^{+}}}{\longrightarrow}(\xi-1)^{-\mu_{+}}\mbox{exp}\left(\frac{\lambda_{+}}{\xi-1}\right)$
(74)

with $\mu_{+}$ and $\lambda_{+}$ defined in (73).

  2. 2.

Asymptotic behaviour as $\xi\to+\infty$: By means of the transformation
$\eta=1/\xi$, it is not difficult to verify that the point at infinity is
again an irregular singular point of rank one. Therefore, the asymptotic
behaviour of the solutions to equation (66) can be derived according to the
method outlined in Olver (1994). To this purpose, we observe that

$\mathfrak{p}(\xi)=\sum_{\kappa=0}^{\infty}\frac{\widehat{\mathfrak{f}}_{\kappa}}{\xi^{k}}=\mathcal{O}\left(\frac{1}{\xi^{2}}\right),\quad\mathfrak{q}(\xi)=\sum_{\kappa=0}^{\infty}\frac{\widehat{\mathfrak{g}}_{\kappa}}{\xi^{k}}=4x_{e}^{2}\Omega^{2}+\frac{8x_{e}\Omega^{2}}{\xi}+\mathcal{O}\left(\frac{1}{\xi^{2}}\right).$
(75)

With the help of (26), we immediately find that $\lambda_{\pm}=\pm
2ix_{e}\Omega$ and $\mu_{\pm}=\pm 2i\Omega$. Hence, the QNM boundary condition
at space-like infinity can be expressed as

$\psi_{\Omega\ell\epsilon}\underset{{\xi\to+\infty}}{\longrightarrow}\xi^{2i\Omega}e^{2ix_{e}\Omega\xi}.$
(76)

Let us transform the radial function $\psi_{\Omega\ell\epsilon}(\xi)$ into a
new radial function $\Phi_{\Omega\ell\epsilon}(\xi)$ such that the QNM
boundary conditions are automatically implemented and $\Phi_{\sigma\ell
s}(\xi)$ is regular at $\xi=1$ and at space-like infinity. To this aim, we
consider the transformation

$\psi_{\Omega\ell\epsilon}(\xi)=\xi^{2i\Omega+\mu_{+}}(\xi-1)^{-\mu_{+}}e^{2ix_{e}\Omega(\xi-1)+\frac{\lambda_{+}}{\xi-1}}\Phi_{\Omega\ell\epsilon}(\xi).$
(77)

If we rewrite it in a more compact form, namely

$\displaystyle\psi_{\Omega\ell\epsilon}(\xi)$ $\displaystyle=$
$\displaystyle\xi^{2ia\Omega}(\xi-1)^{-2i(a-1)\Omega}e^{2ix_{e}\Omega\eta(\xi)}\Phi_{\Omega\ell\epsilon}(\xi),$
(78) $\displaystyle\eta(\xi)$ $\displaystyle=$
$\displaystyle\xi-1+\frac{1}{(\mu_{e}^{2}x_{e}^{2}-1)(\xi-1)},\quad
a=1+\frac{2x_{e}(\mu_{e}^{4}x_{e}^{4}-\mu_{e}^{2}x_{e}^{2}-1)}{3(\mu_{e}^{2}x_{e}^{2}-1)^{2}},$
(79)

and we replace it into (64), we end up with the differential equation

$P_{2e}(\xi)\Phi^{{}^{\prime\prime}}_{\Omega\ell\epsilon}(\xi)+P_{1e}(\xi)\Phi^{{}^{\prime}}_{\Omega\ell\epsilon}(\xi)+P_{0e}(\xi)\Phi_{\Omega\ell\epsilon}(\xi)=0$
(80)

with

$\displaystyle P_{2e}(\xi)$ $\displaystyle=$
$\displaystyle\frac{\xi^{2}(\xi-1)^{2}}{4}f^{2}_{e}(\xi),$ (81) $\displaystyle
P_{1e}(\xi)$ $\displaystyle=$
$\displaystyle\xi(\xi-1)f_{e}(\xi)\left\\{\frac{\xi(\xi-1)}{4}f^{{}^{\prime}}_{e}(\xi)+i\Omega
f_{e}(\xi)\left[x_{e}\xi(\xi-1)\eta^{{}^{\prime}}(\xi)+\xi-a\right]\right\\},$
(82) $\displaystyle P_{0e}(\xi)$ $\displaystyle=$
$\displaystyle-\mathfrak{Q}_{+}(\xi)\mathfrak{Q}_{-}(\xi)\Omega^{2}+i\Omega
f_{e}(\xi)\mathfrak{L}(\xi)-\xi^{2}(\xi-1)^{2}V_{\epsilon}(\xi),$ (83)
$\displaystyle\mathfrak{Q}_{\pm}(\xi)$ $\displaystyle=$ $\displaystyle
x_{e}\xi(\xi-1)f_{e}(\xi)\eta^{{}^{\prime}}(\xi)\pm
x_{e}\xi(\xi-1)+(\xi-a)f_{e}(\xi),$ (84) $\displaystyle\mathfrak{L}(\xi)$
$\displaystyle=$
$\displaystyle\frac{x_{e}}{2}\xi^{2}(\xi-1)^{2}\left[f_{e}(\xi)\eta^{{}^{\prime}}(\xi)\right]^{{}^{\prime}}+\frac{1}{2}\xi(\xi-1)(\xi-a)f_{e}^{{}^{\prime}}(\xi)-\frac{1}{2}(\xi^{2}-2a\xi+a)f_{e}(\xi).$
(85)

Let us introduce the transformation $\xi=2/(1-y)$ mapping the point at
infinity and the event horizon to $y=1$ and $y=-1$, respectively. Furthermore,
a dot denotes differentiation with respect to the new variable $y$. Then,
equation (80) becomes

$S_{2e}(y)\ddot{\Phi}_{\Omega\ell\epsilon}(y)+S_{1e}(y)\dot{\Phi}_{\Omega\ell\epsilon}(y)+S_{0e}(y)\Phi_{\Omega\ell\epsilon}(y)=0,$
(86)

where

$\displaystyle S_{2e}(y)$ $\displaystyle=$
$\displaystyle\frac{(1+y)^{2}}{4}f_{e}^{2}(y),$ (87) $\displaystyle S_{1e}(y)$
$\displaystyle=$ $\displaystyle
i\Omega\frac{1+y}{1-y}f_{e}^{2}(y)\left[x_{e}(1-y^{2})\dot{\eta}(y)+2-a(1-y)\right]-\frac{(1+y)^{2}}{2(1-y)}f_{e}^{2}(y)+\frac{(1+y)^{2}}{4}f_{e}(y)\dot{f}_{e}(y),$
(88) $\displaystyle S_{0e}(y)$ $\displaystyle=$
$\displaystyle\Omega^{2}\Sigma_{2e}(y)+i\Omega\Sigma_{1e}(y)+\Sigma_{0e}(y)$
(89)

with

$\displaystyle\Sigma_{2e}(y)$ $\displaystyle=$
$\displaystyle\frac{4x_{e}^{2}(1+y)^{2}}{(1-y)^{4}}-\frac{f_{e}^{2}(y)}{(1-y)^{2}}\left[x_{e}(1-y^{2})\dot{\eta}(y)+2-a(1-y)\right]^{2},$
(90) $\displaystyle\Sigma_{1e}(y)$ $\displaystyle=$
$\displaystyle\frac{x_{e}}{2}(1+y)^{2}f_{e}(y)\dot{f}_{e}(y)\dot{\eta}(y)+\frac{x_{e}}{2}\frac{(1+y)^{2}}{1-y}f_{e}^{2}(y)\left[(1-y)\ddot{\eta}(y)-2\dot{\eta}(y)\right]+$
(91)
$\displaystyle\frac{1+y}{2(1-y)}\left[2-a(1-y)\right]f_{e}(y)\dot{f}_{e}(y)-\frac{f_{e}^{2}(y)}{(1-y)^{2}}\left[2-2a(1-y)+\frac{a}{2}(1-y)^{2}\right],$
$\displaystyle\Sigma_{1e}(y)$ $\displaystyle=$
$\displaystyle-\frac{4(1+y)^{2}}{(1-y)^{4}}V_{\epsilon}(y).$ (92)

and the requirement that $\Phi_{\Omega\ell\epsilon}(y)$ is regular at $y=\pm
1$. As a result of the transformation introduced above, we have

$f_{e}(y)=1-\frac{1-y}{2x_{e}}\mbox{erf}\left(\frac{2\mu_{e}x_{e}}{1-y}\right)+\frac{2\mu_{e}}{\sqrt{\pi}}e^{-\frac{4\mu_{e}^{2}x_{e}^{2}}{(1-y)^{2}}},\qquad
V_{\epsilon}(y)=\frac{(1-y)^{2}}{16}f_{e}(y)\left[\epsilon(1-y)\dot{f}_{e}(y)+\ell(\ell+1)\right]$
(93)

and

$\eta(y)=\frac{1+y}{1-y}+\frac{1-y}{(\mu_{e}^{2}x_{e}^{2}-1)(1+y)}.$ (94)

Table 3: Classification of the points $y=\pm 1$ for the relevant functions entering in (87), (88) and (89). The abbreviations $z$ ord $n$ and $p$ ord $m$ stand for zero of order $n$ and pole of order $m$, respectively. $y$ | $f_{e}(y)$ | $V_{\epsilon}(y)$ | $\eta(y)$ | $S_{2e}(y)$ | $S_{1e}(y)$ | $S_{0e}(y)$
---|---|---|---|---|---|---
$-1$ | z ord 2 | z ord 2 | p ord 1 | z ord 6 | z ord 4 | z ord 4
$+1$ | $+1$ | z ord 2 | p ord 1 | $+1$ | p ord 2 | p ord 2

Table 3 shows that the coefficients of the differential equation (86) share a
common zero of order $4$ at $y=-1$ while $y=1$ is a pole of order $2$ for the
coefficients $S_{1e}(y)$ and $S_{0e}(y)$. Hence, in order to apply the
spectral method, we need to multiply (86) by $(1-y)^{2}/(1+y)^{4}$. As a
result, we end up with the following differential equation

$M_{2e}(y)\ddot{\Phi}_{\Omega\ell\epsilon}(y)+M_{1e}(y)\dot{\Phi}_{\Omega\ell\epsilon}(y)+M_{0e}(y)\Phi_{\Omega\ell\epsilon}(y)=0,$
(95)

where

$M_{2e}(y)=\frac{(1-y)^{2}}{4(1+y)^{2}}f_{e}^{2}(y),\quad M_{1e}(y)=i\Omega
N_{1e}(y)+N_{0e}(y),\quad M_{0e}(y)=\Omega^{2}C_{2e}(y)+i\Omega
C_{1e}(y)+C_{0e}(y)$ (96)

with

$\displaystyle N_{1e}(y)$ $\displaystyle=$
$\displaystyle\frac{1-y}{(1+y)^{3}}f_{e}^{2}(y)\left[x_{e}(1-y^{2})\dot{\eta}(y)+2-a(1-y)\right],\quad
N_{0e}(y)=\frac{f_{e}(y)}{4(1+y)^{2}}\frac{d}{dy}\left((1-y)^{2}f_{e}(y)\right),$
(97) $\displaystyle C_{2e}(y)$ $\displaystyle=$
$\displaystyle\frac{4x_{e}^{2}}{(1-y^{2})^{2}}-\frac{f_{e}^{2}(y)}{(1+y)^{4}}\left[x_{e}(1-y^{2})\dot{\eta}(y)+2-a(1-y)\right]^{2},$
(98) $\displaystyle C_{1e}(y)$ $\displaystyle=$
$\displaystyle\frac{x_{e}}{2}\left(\frac{1-y}{1+y}\right)^{2}f_{e}(y)\dot{f}_{e}(y)\dot{\eta}(y)+\frac{x_{e}}{2}\frac{1-y}{(1+y)^{2}}f_{e}^{2}(y)\left[(1-y)\ddot{\eta}(y)-2\dot{\eta}(y)\right]+$
(99)
$\displaystyle\frac{1-y}{2(1+y)^{3}}\left[2-a(1-y)\right]f_{e}(y)\dot{f}_{e}(y)-\frac{f_{e}^{2}(y)}{(1+y)^{4}}\left[2-2a(1-y)+\frac{a}{2}(1-y)^{2}\right],$
$\displaystyle C_{0e}(y)$ $\displaystyle=$
$\displaystyle-\frac{4V_{\epsilon}(y)}{(1-y^{2})^{2}}.$ (100)

It can be easily checked with Maple that

$\displaystyle\lim_{y\to
1^{-}}M_{2e}(y)=0=\lim_{y\to-1^{+}}M_{2e}(y),\quad\lim_{y\to
1^{-}}M_{1e}(y)=\frac{1}{2}ix_{e}\Omega,\quad\lim_{y\to-1^{+}}M_{1e}(y)=i\Omega\Lambda_{1e},$
(101) $\displaystyle\lim_{y\to
1^{-}}M_{0e}(y)=A_{2e}\Omega^{2}+A_{0e},\quad\lim_{y\to-1^{+}}M_{0e}(y)=B_{2e}\Omega^{2}+B_{0e},$
(102)

where

$\displaystyle\Lambda_{1e}$ $\displaystyle=$
$\displaystyle-\frac{x_{e}}{2}\left(\mu_{e}^{2}x_{e}^{2}-1\right),\quad
A_{2e}=\frac{2\mu_{e}^{4}x_{e}^{6}+3\mu_{e}^{4}x_{e}^{4}+\mu_{e}^{2}x_{e}^{4}-6\mu_{e}^{2}x_{e}^{2}-5x_{e}^{2}+3}{6\left(\mu_{e}^{2}x_{e}^{2}-1\right)^{2}},\quad
A_{0e}=-\frac{\ell(\ell+1)}{16},$ (103) $\displaystyle B_{2e}$
$\displaystyle=$
$\displaystyle\frac{x_{e}\left(2\mu_{e}^{8}x_{e}^{9}+35\mu_{e}^{6}x_{e}^{7}+18\mu_{e}^{6}x_{e}^{6}-107\mu_{e}^{4}x_{e}^{5}-54\mu_{e}^{4}x_{e}^{4}+88\mu_{e}^{2}x_{e}^{3}+54\mu_{e}^{2}x_{e}^{2}-10x_{e}-18\right)}{36\left(\mu_{e}^{2}x_{e}^{2}-1\right)^{2}},$
(104) $\displaystyle B_{0e}$ $\displaystyle=$
$\displaystyle-\frac{\mu_{e}^{2}x_{e}^{2}-1}{16}\ell(\ell+1).$ (105)

Finally, in order to apply the spectral method, we rewrite the differential
equation (95) into the following form

$\widehat{L}^{(e)}_{0}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]+i\widehat{L}^{(e)}_{1}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]\Omega+\widehat{L}_{2}^{(e)}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]\Omega^{2}=0$
(106)

with

$\displaystyle\widehat{L}^{(e)}_{0}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]$
$\displaystyle=$
$\displaystyle\widehat{L}^{(e)}_{00}(y)\Phi_{\Omega\ell\epsilon}+\widehat{L}^{(e)}_{01}(y)\dot{\Phi}_{\Omega\ell\epsilon}+\widehat{L}^{(e)}_{02}(y)\ddot{\Phi}_{\Omega\ell\epsilon},$
(107)
$\displaystyle\widehat{L}^{(e)}_{1}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]$
$\displaystyle=$
$\displaystyle\widehat{L}^{(e)}_{10}(y)\Phi_{\Omega\ell\epsilon}+\widehat{L}^{(e)}_{11}(y)\dot{\Phi}_{\Omega\ell\epsilon}+\widehat{L}^{(e)}_{12}(y)\ddot{\Phi}_{\Omega\ell\epsilon},$
(108)
$\displaystyle\widehat{L}^{(e)}_{2}\left[\Phi_{\Omega\ell\epsilon},\dot{\Phi}_{\Omega\ell\epsilon},\ddot{\Phi}_{\Omega\ell\epsilon}\right]$
$\displaystyle=$
$\displaystyle\widehat{L}^{(e)}_{20}(y)\Phi_{\Omega\ell\epsilon}+\widehat{L}^{(e)}_{21}(y)\dot{\Phi}_{\Omega\ell\epsilon}+\widehat{L}^{(e)}_{22}(y)\ddot{\Phi}_{\Omega\ell\epsilon}.$
(109)

Moreover, in Table 4, we have summarized the $\widehat{L}^{(e)}_{ij}$
appearing in (107)–(109) and their limiting values at $y=\pm 1$.

Table 4: Definitions of the coefficients $\widehat{L}^{(e)}_{ij}$ and their corresponding behaviours at the endpoints of the interval $-1\leq y\leq 1$. The symbols appearing in this table have been defined in (96)-(105). $(i,j)$ | $\displaystyle{\lim_{y\to-1^{+}}}\widehat{L}^{(e)}_{ij}$ | $\widehat{L}^{(e)}_{ij}$ | $\displaystyle{\lim_{y\to 1^{-}}}\widehat{L}^{(e)}_{ij}$
---|---|---|---
$(0,0)$ | $B_{0e}$ | $C_{0e}$ | $A_{0e}$
$(0,1)$ | $0$ | $N_{0e}$ | $0$
$(0,2)$ | $0$ | $M_{2e}$ | $0$
$(1,0)$ | $0$ | $C_{1e}$ | $0$
$(1,1)$ | $\Lambda_{1e}$ | $N_{1e}$ | $x_{e}/2$
$(1,2)$ | $0$ | $0$ | $0$
$(2,0)$ | $B_{2e}$ | $C_{2e}$ | $A_{2e}$
$(2,1)$ | $0$ | $0$ | $0$
$(2,2)$ | $0$ | $0$ | $0$

##  III Numerical method

In order to solve the differential eigenvalue problem (106) to determine the
QNMs along with the corresponding frequencies $\Omega$, we have to discretize
the differential operators $\widehat{L}^{(e)}_{j}[\cdot]$ with $j=1,2,3$
defined in (107)-(109). Since our problem is posed on the finite interval
$[-1,1]$ without any boundary conditions, more precisely, we only require that
the QNMs be regular functions at $y=\pm 1$, then, it is natural to choose a
Tchebyshev-type spectral method Trefethen (2000); Boyd (2000). Namely, we are
going to expand the function $y\mapsto\Phi_{\Omega\ell\epsilon}(y)$ in the
form of a truncated Tchebyshev series

$\Phi_{\Omega\ell\epsilon}(y)=\sum_{k=0}^{N}a_{k}T_{k}(y),$ (110)

where $N\ \in\ \mathbb{N}$ is kept as a numerical parameter,
$\\{a_{k}\\}_{k=0}^{N}\ \subseteq\ \mathds{R}$, and $\\{T_{k}(y)\\}_{k=0}^{N}$
are the Tchebyshev polynomials of the first kind

$T_{k}:[-1,1]\ \longrightarrow\ [-1,1]\,,\qquad y\ \longmapsto\
\cos\,\bigl{(}k\arccos y\bigr{)}\,.$ (111)

After substituting expansion (110) into the differential equation (106), we
obtain an eigenvalue problem with polynomial coefficients. In order to
translate it into the realm of numerical linear algebra, we employ the
collocation method Boyd (2000). Specifically, rather than insisting that the
polynomial function in $y$ is identically zero (a condition equivalent to
having polynomial solutions for the differential problem as per equation
(106)), we impose a weaker requirement. This involves ensuring that the
polynomial vanishes at $N+1$ strategically selected points. The number $N+1$
coincides exactly with the number of unknown coefficients
$\\{a_{k}\\}_{k=0}^{N}$. For the collocation points, we implemented the
Tchebyshev roots grid Fox and Parker (1968)

$y_{k}=-\cos{\left(\frac{(2k+1)\pi}{2(n+1)}\right)},\quad
k\in\\{0,1,\ldots,N\\}.$ (112)

In our numerical codes, we also implemented the second option of the
Tchebyshev extrema grid

$y_{k}=-\cos{\left(\frac{k\pi}{n}\right)},\quad k\in\\{0,1,\ldots,N\\}.$ (113)

The users are free to choose their favourite collocation points. Notice that
we used the roots grid in our computation, and in any case, the theoretical
performance of the two available options is known to be absolutely comparable
Fox and Parker (1968); Boyd (2000).

Upon implementing the collocation method, we derive a classical matrix-based
quadratic eigenvalue problem, as detailed in Tisseur and Meerbergen (2001)

$(M_{0}+iM_{1}\Omega+M_{2}\Omega^{2})\bf{a}=\bf{0}.$ (114)

In this formulation, the square real matrices $M_{j}$, each of size
$(N+1)\times(N+1)$ for $j=0,1,2$, represent the spectral discretizations of
the operators $\widehat{L}^{(e)}_{j}[\cdot]$, respectively. The problem (114)
is solved numerically with the polyeig function from Matlab. This polynomial
eigenvalue problem yields $2(N+1)$ potential values for the parameter
$\Omega$. To discern the physical values of $\Omega$ that correspond to the
black hole’s QNMs, we first overlap the root plots for various values of $N$
in equation (110), such as $N=100,150,200$. We then identify the consistent
roots whose positions remain stable across these different $N$ values.

In order to reduce the rounding and other floating point errors, we performed
all our computations with multiple precision arithmetic that is built in Maple
and which is brought into Matlab by the Advanpix toolbox Holodoborodko (2023).
All numerical computations reported in this study have been performed with
$200$ decimal digits accuracy.

##  IV Numerical results

In this section, we delve into the analysis of the QNMs for both non-extreme
and extreme noncommutative geometry-inspired Schwarzschild black holes. The
central revelation from our numerical computations is the stability of both
black hole regimes under scalar, electromagnetic, and gravitational
perturbations.

The fact that the noncommutative geometry-inspired Schwarzschild black hole
reproduces its classical counterpart for large values of the rescaled mass
parameter $\mu$ provides a basis for validating our numerical method. By
selecting a sufficiently large $\mu$-value, computing the QNMs, and
demonstrating their agreement with the results obtained by Leaver (1985); Iyer
(1987); Mamani et al. (2022) for the classic Schwarzschild black hole in the
$s=0$ case, we can substantiate the accuracy of our approach. To this end, we
refer to Table 5, where the numerical values listed in the last column were
derived using our high-precision computational approach, employing 200
polynomials with an accuracy of $200$ digits.

As an additional note, we would like to highlight that in the extreme case,
the transition to the Schwarzschild limit is more complex due to potential
non-analyticities in the behavior of the metric and associated physical
quantities as $\theta$ tends towards zero. The extreme configuration in
noncommutative geometry, characterized by a specific critical mass $\mu_{e}$
at which the inner and outer horizons converge, does not straightforwardly
reduce to the Schwarzschild solution by merely setting $\theta$ to zero. This
non-analyticity arises because the extreme condition itself alters the
underlying spacetime structure in a way that does not directly map onto the
classical Schwarzschild case without noncommutativity.

Table 5: This table presents the quasinormal frequencies for scalar perturbations ($s=0$) of the classic Schwarzschild black hole alongside those for a noncommutative geometry-inspired Schwarzschild black hole with a large mass parameter (refer to the last column). The third column showcases numerical values determined by Leaver (1985) using the continued fraction method. The fourth column includes third-order WKB results reported by Iyer (1987), and the fifth column presents numerical values from Mamani et al. (2022), derived via the spectral method employing 40 polynomials. Here, $\Omega$ represents the dimensionless frequency as defined in equation (10). The notation ’N/A’ indicates data not available. $\ell$ | $n$ |  $\Omega$ Leaver (1985) |  $\Omega$ Iyer (1987) |  $\Omega$ Mamani et al. (2022) |  $\Omega$, $\mu=10^{3}$
---|---|---|---|---|---
$0$ | $0$ | $0.1105-0.1049i$ | $0.1046-0.1152i$ | $0.1105-0.1049i$ | $0.1104549-0.1048957i$
| $1$ | $0.0861-0.3481i$ | $0.0892-0.3550i$ | N/A | $0.0861169-0.3480525i$
$1$ | $0$ | $0.2929-0.0977i$ | $0.2911-0.0980i$ | $0.2929-0.0977i$ | $0.2929361-0.0976599i$
| $1$ | $0.2645-0.3063i$ | $0.2622-0.3704i$ | $0.2645-0.3063i$ | $0.2644487-0.3062574i$
| $2$ | $0.2295-0.5401i$ | $0.2235-0.5268i$ | N/A | $0.2295393-0.5401334i$
| $3$ | $0.2033-0.7883i$ | $0.1737-0.7486i$ | N/A | $0.2032584-0.7882978i$
$2$ | $0$ | $0.4836-0.0968i$ | $0.4832-0.0968i$ | $0.4836-0.0968i$ | $0.4836439-0.0967588i$
| $1$ | $0.4639-0.2956i$ | $0.4632-0.2958i$ | $0.4639-0.2956i$ | $0.4638506-0.2956039i$
| $2$ | $0.4305-0.5086i$ | $0.4317-0.5034i$ | $0.4305-0.5086i$ | $0.4305441-0.5085584i$
| $3$ | $0.3939-0.7381i$ | $0.3926-0.7159i$ | N/A | $0.3938631-0.7380966i$

It is interesting to observe that Batic et al. (2019b) utilized a sixth-order
WKB approximation to show that the WKB method does not converge in critical
cases marked by third-order instabilities. These occur for specific choices of
the rescaled mass parameters close to $\mu=\mu_{e}$ for a nonextreme,
noncommutative geometry-inspired Schwarzschild black hole. Instabilities were
also detected by Batic et al. (2019b) in the case of an extreme black hole.
This indicates that conventional methods have exhibited inconsistencies,
particularly when applied to non-commutative geometry settings. There are
several reasons for this. For instance, the WKB approximation relies on the
local properties of the potential near its peak. However, it fails in cases
where the second derivative of the potential at the peak is zero. As for the
Inverse Potential and Asymptotic Iteration Methods, the former depends on a
rough approximation of the potential at the maximum, while the latter depends
on the choice of initial parameters or specific local conditions that may not
be optimally representative or easily determinable in complex geometries like
those modified by the non-commutative parameter. These dependencies can lead
to instabilities or divergences in computations, as indicated by the lack of
convergence observed in higher-order WKB calculations, suggesting that these
methods may reach their analytical limitations in such scenarios.

Unlike local methods, the spectral method utilizes the global properties of
the potential, thus avoiding inaccuracies due to local approximations.
Moreover, it provides high accuracy and exponential convergence rates for the
solutions of differential equations once the irregular part of the solution
has been analytically captured. This is particularly beneficial in the quest
for QNMs of astrophysical objects. Furthermore, the method is less sensitive
to the choice of initial parameters (in contrast with the Asymptotic Iterative
Method), making it more robust and reliable for theoretical explorations in
gravitational physics where exact solutions are often challenging to compute.
For non-commutative geometry-inspired black holes, the spectral method allows
us to capture the modified dynamics without the constraints and approximations
imposed by traditional techniques. This capability is crucial in revealing
that previously reported instabilities are likely artifacts of the limitations
of those methods rather than inherent physical instabilities in the black hole
models themselves.

Consequently, by harnessing the global analytical capabilities of the spectral
method, we have circumvented the shortcomings of traditional approaches like
the WKB method. This has allowed us to demonstrate that the instabilities
previously reported in non-commutative geometry-inspired black holes are not
inherent to the black holes themselves but are rather artifacts resulting from
the limitations of local approximation methods.

Therefore, our application of the spectral method provides a clearer, more
accurate depiction of the stability characteristics of these black holes. For
a comprehensive presentation of our numerical results on massless scalar,
electromagnetic, and gravitational perturbations, we invite readers to consult
Tables 6, 7, and 8. These tables detail the calculated QNMs of a nonextreme
noncommutative black hole, illustrating variations across different values of
the mass parameter. In contrast, Tables 9, 10, and 11 focus on the extreme
case, providing analogous data for the same spectrum of perturbations.

Table 6: QNMs for scalar perturbations ($s=0$) of the noncommutative geometry inspired Schwarzschild metric (nonextreme case) for different values of the mass parameter $\mu$. The third and fifth columns represent the numerical values found by Batic et al. (2019b) by means of a sixth-order WKB approximation. The fourth and sixth columns report the corresponding results obtained by means of the spectral method employing $150$ polynomials with an accuracy of $200$ digits. Here, $\Omega$ represents the dimensionless frequency as defined in equation (10). The ’SM’ abbreviation denotes the ’Spectral Method’. $\ell$ | $n$ |  $\Omega$, $\mu=1.95$ Batic et al. (2019b) |  $\Omega$, $\mu=1.95$ (SM) |  $\Omega$, $\mu=2.25$ Batic et al. (2019b) |  $\Omega$, $\mu=2.25$ (SM)
---|---|---|---|---|---
$0$ | $0$ | $0.0397+0.1330i$ | $0.1041-0.0930i$ | $0.0222-0.0944i$ | $0.1101-0.1010i$
| $1$ | $0.1555+0.4548i$ | $0.0343-0.3374i$ | $0.0880+0.3947i$ | $0.0595-0.3298i$
$1$ | $0$ | $0.2727-0.0818i$ | $0.2890-0.0908i$ | $0.2882-0.3509i$ | $0.2916-0.0958i$
| $1$ | $0.1256-0.2689i$ | $0.2381-0.2807i$ | $0.2449-0.3059i$ | $0.2556-0.2971i$
| $2$ | $0.0779+0.5589i$ | $0.1480-0.5047i$ | $0.1771-0.5284i$ | $0.1961-0.5186i$
| $3$ | $0.3072+0.8805i$ | $0.0611-0.7768i$ | $0.0862-0.7594i$ | $0.1056-0.7762i$
$2$ | $0$ | $0.4767-0.0889i$ | $0.4806-0.0916i$ | $0.4819-0.0964i$ | $0.4825-0.0957i$
| $1$ | $0.4209-0.2669i$ | $0.4477-0.2667i$ | $0.4561-0.2944i$ | $0.4581-0.2909i$
| $2$ | $0.2999-0.4590i$ | $0.3833-0.4607i$ | $0.4122-0.5015i$ | $0.4134-0.4965i$
| $3$ | $0.1285-0.6908i$ | $0.2925-0.6865i$ | $0.3544-0.7154i$ | $0.3537-0.7160i$

Our findings suggest a remarkable stability of the nonextreme noncommutative
Schwarzschild black hole across various perturbations. This observation is
particularly noteworthy given the existing literature, such as the study by
Brown and Mann (2011), which posited potential instability in similar
scenarios. Our results, however, indicate that even in extreme scenarios where
the Cauchy horizon coincides with the event horizon, the black hole exhibits
stability under the class of perturbations considered. This result extends the
work of Batic and Nicolini (2010), where the stability of the nonextreme
noncommutative Schwarzschild black hole’s interior was established by
investigating the propagation of a massless scalar field between the two
horizons.

Table 7: QNMs for electromagnetic perturbations ($s=1$) of the noncommutative geometry inspired Schwarzschild metric (nonextreme case) for different values of the mass parameter $\mu$. The corresponding numerical results are obtained by means of the spectral method employing $150$ polynomials with an accuracy of $200$ digits. Here, $\Omega$ represents the dimensionless frequency as defined in equation (10). $\ell$ | $n$ |  $\Omega$, $\mu=1.95$ |  $\Omega$, $\mu=2.25$
---|---|---|---
$1$ | $0$ | $0.2469-0.0831i$ | $0.2478-0.0899i$
| $1$ | $0.1945-0.2569i$ | $0.2090-0.2810i$
| $2$ | $0.0965-0.4661i$ | $0.1491-0.4936i$
| $3$ | $0.0181-0.6972i$ | $0.0656-0.7133i$
$2$ | $0$ | $0.4552-0.0887i$ | $0.4566-0.0936i$
| $1$ | $0.4220-0.2682i$ | $0.4314-0.2848i$
| $2$ | $0.3566-0.4561i$ | $0.3855-0.4866i$
| $3$ | $0.2637-0.6656i$ | $0.3254-0.7022i$
Table 8: QNMs for vector-type gravitational perturbations ($s=2$) of the noncommutative geometry inspired Schwarzschild metric (nonextreme case) for different values of the mass parameter $\mu$. The corresponding numerical results are obtained by means of the spectral method employing $150$ polynomials with an accuracy of $200$ digits. Here, $\Omega$ represents the dimensionless frequency as defined in equation (10). $\ell$ | $n$ |  $\Omega$, $\mu=1.95$ |  $\Omega$, $\mu=2.25$
---|---|---|---
$2$ | $0$ | $0.3723-0.0776i$ | $0.3726-0.0861i$
| $1$ | $0.3411-0.2316i$ | $0.3437-0.2611i$
| $2$ | $0.2862-0.3945i$ | $0.2973-0.4462i$
| $3$ | $0.2073-0.5902i$ | $0.2477-0.6492i$
$3$ | $0$ | $0.5979-0.0857i$ | $0.5985-0.0913i$
| $1$ | $0.5730-0.2570i$ | $0.5786-0.2755i$
| $2$ | $0.5252-0.4298i$ | $0.5417-0.4649i$
| $3$ | $0.4569-0.6104i$ | $0.4923-0.6632i$
| $4$ | $0.3693-0.8090i$ | $0.4347-0.8720i$
Table 9: This table presents the quasinormal frequencies for scalar perturbations ($s=0$) of the extreme noncommutative geometry-inspired Schwarzschild black hole. The third and fourth columns showcase numerical values determined by Batic et al. (2019b) using the WKB approximation and the inverted potential method (PT), respectively. The fifth column includes the results obtained via the spectral method employing $150$ polynomials with an accuracy of $200$ digits. Here, $\Omega$ represents the dimensionless frequency as defined in equation (64). The ’SM’ abbreviation denotes the ’Spectral Method’. $\ell$ | $n$ |  $\Omega$ (WKB) Batic et al. (2019b) |  $\Omega$ (PT) Batic et al. (2019b) |  $\Omega$ (SM)
---|---|---|---|---
$0$ | $0$ | $0.0395+0.1367i$ | $0.1003-0.1232i$ | $0.1027-0.0928i$
| $1$ | $0.1524+0.4582i$ | $0.1003-0.3697i$ | $0.0299-0.3374i$
$1$ | $0$ | $0.2699-0.0744i$ | $0.2983-0.0997i$ | $0.2882-0.0895i$
| $1$ | $0.1061-0.2684i$ | $0.2983-0.2991i$ | $0.2341-0.2777i$
| $2$ | $0.1023+0.5721i$ | $0.2983-0.4985i$ | $0.1400-0.5021i$
| $3$ | $0.3374+0.9027i$ | $0.2983-0.6979i$ | $0.0317-0.5720i$
$2$ | $0$ | $0.4756-0.0870i$ | $0.4881-0.0949i$ | $0.4802-0.0903i$
| $1$ | $0.4136-0.2612i$ | $0.4881-0.2847i$ | $0.4447-0.2731i$
| $2$ | $0.2782-0.4540i$ | $0.4881-0.4745i$ | $0.3760-0.4657i$
| $3$ | $0.0918-0.6961i$ | $0.4881-0.6643i$ | $0.2806-0.6818i$
Table 10: This table presents the quasinormal frequencies for electromagnetic perturbations ($s=1$) of the extreme noncommutative geometry-inspired Schwarzschild black hole via the spectral method employing $150$ polynomials with an accuracy of $200$ digits. Here, $\Omega$ represents the dimensionless frequency as defined in equation (64). $\ell$ | $n$ | $\Omega$
---|---|---
$1$ | $0$ | $0.2465-0.0813i$
| $1$ | $0.1904-0.2527i$
| $2$ | $0.0883-0.4627i$
$2$ | $0$ | $0.4549-0.0872i$
| $1$ | $0.4192-0.2640i$
| $2$ | $0.3493-0.4503i$
| $3$ | $0.2516-0.6601i$
Table 11: This table presents the quasinormal frequencies for vector-type gravitational perturbations ($s=2$) of the extreme noncommutative geometry-inspired Schwarzschild black hole via the spectral method employing $150$ polynomials with an accuracy of $200$ digits. Here, $\Omega$ represents the dimensionless frequency as defined in equation (64). $\ell$ | $n$ | $\Omega$
---|---|---
$2$ | $0$ | $0.3723-0.0749i$
| $1$ | $0.3394-0.2240i$
| $2$ | $0.2811-0.3853i$
| $3$ | $0.1975-0.5832i$
$3$ | $0$ | $0.5979-0.0839i$
| $1$ | $0.5711-0.2516i$
| $2$ | $0.5199-0.4214i$
| $3$ | $0.4471-0.6005i$
| $4$ | $0.3546-0.7997i$

##  V Conclusions and outlook

In this study, we have presented a comprehensive analysis of QNMs for
noncommutative geometry-inspired Schwarzschild black holes, encompassing both
non-extreme and extreme cases. Our investigation, employing the spectral
method, covered scalar, electromagnetic, and gravitational perturbations. The
significance of our findings is the remarkable stability exhibited by both
non-extreme and extreme versions of these black holes across various
perturbations. Contrary to previous claims in the literature, our results
suggest that the noncommutative geometry-inspired Schwarzschild black holes
are stable, a finding that is crucial for the understanding of their physical
properties. This is particularly relevant in light of the existing debate over
the stability of such black holes. Our study, therefore, provides a new
perspective on the dynamics of these exotic objects and contributes to a more
nuanced understanding of their behaviour under perturbations.

###  V.1 Perspectives

Our future endeavours will focus on expanding this study in several key
directions. Firstly, we plan to investigate the QNMs in the context of massive
perturbations. This extension will allow us to explore the influence of mass
on the stability and dynamics of these black holes, offering further insights
into their complex nature. Additionally, we aim to apply our methods to the
study of QNMs of noncommutative geometry-inspired wormholes Garattini and Lobo
(2009); Nicolini and Spallucci (2010). This will not only broaden our
understanding of wormhole physics but also provide a comparative framework to
evaluate the similarities and differences between black holes and wormholes
within the context of noncommutative geometry. Another intriguing direction
for future research involves the investigation of the so-called dirty black
holes as proposed by Nicolini and Spallucci (2010). These black holes,
characterized by additional matter fields surrounding them, present a more
realistic and complex scenario. Analyzing their QNMs will shed light on how
external matter influences the stability and quasinormal spectra of black
holes, further enriching our understanding of these objects. In summary, our
research opens new avenues for exploring the intricate dynamics of
noncommutative geometry-inspired black holes and wormholes. By extending our
analysis to include massive perturbations and dirty black holes, we anticipate
uncovering new aspects of these entities, contributing to the evolving
landscape of gravitational and astrophysical research.

## Code availability

All analytical computations reported in this manuscript have been rechecked in
the computer algebra system Maple. Two Maple sheets corresponding to the
extreme and non-extreme cases can be found in supplementary materials _and_ in
the repository below. The discretization of differential operators (107)–(109)
using the Tchebyshev-type spectral method is equally performed in Maple
computer algebra system. Finally, the numerical solution of the resulting
quadratic eigenvalue problem (114) is performed in Matlab software using the
polyeig function. All these materials are freely accessible at the following
repository

  * •

https://github.com/dutykh/noncommschwarzschild/

## Acknowledgements

This publication is based upon work supported by the Khalifa University of
Science and Technology under Award No. FSU$-2023-014$.

## References

  * Regge and Wheeler (1957) T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957).
  * Hall (2013) B. C. Hall, _Quantum Theory for Mathematicians_ , vol. 267 of _Graduate Texts in Mathematics_ (Springer: New York, 2013).
  * Ferrari and Mashhoon (1984) V. Ferrari and B. Mashhoon, Phys. Rev. D 30, 295 (1984).
  * Leaver (1985) E. W. Leaver, Proc. R. Soc. Lond. A 402, 285 (1985).
  * Leaver (1986) E. W. Leaver, Phys. Rev. D 34, 384 (1986).
  * Nollert and Schmidt (1992) H. P. Nollert and B. G. Schmidt, Phys.Rev. D 45, 2617 (1992).
  * Ciftci et al. (2003) H. Ciftci, R. L. Hall, and N. Saad, J. Phys. A: Math. Gen. 36, 11807 (2003).
  * Cho et al. (2010) H. T. Cho, A. Cornell, J. Doukas, and W. Naylor, Class. Quantum Grav. 27, 155004 (2010).
  * Cho et al. (2012) H. T. Cho, A. Cornell, J. Doukas, T.-R. Huang, and W. Naylor, Adv. Math. Phys. 2012, 281705 (2012).
  * Batic et al. (2018) D. Batic, M. Nowakowski, and K. Redway, Phys. Rev. D 98, 024017 (2018).
  * Macedo (2019) R. P. Macedo, Phys. Rev. D 99, 088501 (2019).
  * Batic et al. (2019a) D. Batic, M. Nowakowski, and K. Redway, Phys. Rev. D 99, 088502 (2019a).
  * Ciftci et al. (2005) H. Ciftci, R. L. Hall, and N. Saad, Phys. Lett. A 340, 388 (2005).
  * Mamani et al. (2022) L. A. H. Mamani, A. D. D. Masa, L. T. Sanches, and V. T. Zanchin, Eur. Phys. J. C 82, 897 (2022).
  * Iyer (1987) H. S. Iyer, Phys. Rev. D 35, 3632 (1987).
  * Daghigh and Green (2012) R. G. Daghigh and M. D. Green, Phys. Rev. D 85, 127501 (2012).
  * Batic and Nowakowski (2023) D. Batic and M. Nowakowski, Math. Methods Appl. Sci. 46, 10833 (2023).
  * Nicolini et al. (2006) P. Nicolini, A. Smailagic, and E. Spallucci, Phys. Lett. B 632, 547 (2006).
  * Bachelot and Motet-Bachelot (1993) A. Bachelot and A. Motet-Bachelot, Ann. Inst. H. Poincare Phys. Theor. 59, 3 (1993).
  * Gundlach et al. (1994) C. Gundlach, R. H. Price, and J. Pullin, Phys. Rev. D 49, 883 (1994).
  * Andersson (1995) N. Andersson, Phys.Rev. D 51, 353 (1995).
  * Andersson (1997) N. Andersson, Phys.Rev. D 55, 468 (1997).
  * Konoplya (2003) R. A. Konoplya, Adv. Theor. Math. Phys. 6, 1135 (2003).
  * Motl and Neitzke (2003) L. Motl and A. Neitzke, Adv. Theor. Math. Phys. 7, 307 (2003).
  * Natário and Schiappa (2004) J. Natário and R. Schiappa, Adv. Theor. Math. Phys. 8, 1001 (2004).
  * Chakrabarti and Gupta (2006) S. K. Chakrabarti and K. S. Gupta, Int. J .Mod. Phys. A 21, 3565 (2006).
  * Giri (2007) P. R. Giri, Int. J. Mod. Phys. A 22, 2047 (2007).
  * Batic et al. (2019b) D. Batic, N. G. Kelkar, M. Nowakowski, and K. Redway, Eur. Phys. J. C 79, 581 (2019b).
  * Liang (2018a) J. Liang, Chinese Phys. Lett. 35, 010401 (2018a).
  * Liang (2018b) J. Liang, Chinese Phys. Lett. 35, 050401 (2018b).
  * Batic and Nicolini (2010) D. Batic and P. Nicolini, Phys.Lett.B 692, 32 (2010).
  * Brown and Mann (2011) E. Brown and R. Mann, Phys. Lett. B 694, 440 (2011).
  * Abramowitz and Stegun (1972) M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions_ (Dover: New York, 1972).
  * Ince (1956) E. L. Ince, _Ordinary Differential Equations_ (Dover: New York, 1956).
  * Bender and Orszag (1999) C. M. Bender and S. A. Orszag, _Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory_ (Springer: New York, 1999).
  * Olver (1994) F. W. J. Olver, Methods Appl. Anal. 1, 1 (1994).
  * Trefethen (2000) L. N. Trefethen, _Spectral methods in MatLab_ (Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2000), URL http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/spectral.html.
  * Boyd (2000) J. P. Boyd, _Chebyshev and Fourier Spectral Methods_ (Dover Publications, New York, 2000), 2nd ed.
  * Fox and Parker (1968) L. Fox and I. B. Parker, _Chebyshev Polynomials in Numerical Analysis_ (Oxford University Press, 1968).
  * Tisseur and Meerbergen (2001) F. Tisseur and K. Meerbergen, SIAM Review 43, 235 (2001), ISSN 0036-1445.
  * Holodoborodko (2023) P. Holodoborodko, _Multiprecision Computing Toolbox for MATLAB 5.2.5.15444_ (Advanpix LLC., Yokohama, Japan, 2023).
  * Garattini and Lobo (2009) R. Garattini and F. S. N. Lobo, Phys. Lett. B 671, 146 (2009).
  * Nicolini and Spallucci (2010) P. Nicolini and E. Spallucci, Class. Quant. Grav. 27, 015010 (2010).