Noncommutative scalar field in the non-extremal Reissner-Nordstr\"om background: QNM spectrum

In our previous work [18] we constructed a model of a noncommutative, charged and massive scalar field based on the angular twist. Then we used this model to analyze the motion of the scalar field in the Reissner-Nordstr\"om black hole background. In particular, we determined the QNM spectrum analytically in the near-extremal limit. To broaden our analysis, in this paper we apply a well defined numerical method, the continued fraction method and calculate the QNM spectrum for a non-extremal Reissner-Nordstr\"om black hole. To check the validity of our analytic calculations, we compare results of the continued fraction method in the near extremal limit with the analytic results obtained in the previous paper. We find that the results are in good agreement. For completeness, we also study the QNM spectrum in the WKB approximation.


Introduction
The study of black hole perturbations has a long history dating back to the work of Regge and Wheeler [1] and Vishveshwara [2]. The initial impetus came from an attempt to analyze the issue of stability of a black hole. Later on, the study was extended to include various types of perturbations in almost all conceivable backgrounds.
With a recent discovery of gravitational waves [3] this line of investigation became ever more important. After being perturbed, black holes return back to their equilibrium by going through a ring down phase, whose most dominant stage is characterized by the long lasting damped oscillations dubbed as quasinormal modes (QNMs) [4]. This phase is characterized by a discrete set of complex frequencies, called QNM frequencies. While the real part of the frequency corresponds to the actual frequency of the wave dynamics, the imaginary part represents the damping factor. For an overview of the subject, we refer the reader to some excellent reviews [5,6,7,8].
Several methods have been devised to determine the quasinormal modes of a black hole. They range from purely analytic methods [9,10] through various semiclassical ones [11,12,13,14] and all down to the purely numerical methods [15].
A way to describe some effects of quantum gravity is to introduce a noncommutative deformation of space-time. The main idea of noncommutative (NC) geometry is that a space-time as we perceive it might possibly be distorted at some relatively high energy scale, which we label as the scale of noncommutativity l N C . The scale of noncommutativity itself may be anything between the TeV scale and the Planck scale. Moreover, if the space-time is indeed modified at the NC scale l N C , then this distortion should somehow be visible in the QNM spectrum of black holes. In other words, with the current capacities as well as with the new possibilities that have opened up by the detection of gravitational waves, a direct comparison of the theoretical results with the observation may signal a presence of a new physics, i.e. may hint toward a presence of an underlying NC space-time structure above some energy threshold.
First results on NC QNM spectrum were published in [16,17]. There, the authors analyze the NC scalar QNM spectrum of the three dimensional BTZ black hole using a direct integration of the equation of motion for the NC scalar field in the BTZ black hole background.
In our recent paper [18] we have investigated a model of NC scalar and gauge fields coupled to a classical background of the Reissner-Nordström (RN) black hole. This model resulted in a master equation governing a behavior of the scalar perturbations in the presence of a noncommutative structure of space-time. The master equation was then analyzed analytically. However, the analytic treatment is only possible for a highly restricted range of parameters corresponding to the near extremal case of the RN black hole. In order to overcome the shortcomings of the analytic treatment, in the present work we analyze the same master equation by a well defined numerical method, the continued fraction method. This method enables us to find solutions for QNM frequencies for a more general set of system parameters. Of course, the results we obtain may also be confronted with the results obtained in [18] after specializing them to the near extremal case. Since the analytic treatment of the reference [18] was carried out under this rather restrictive assumption, the treatment conducted here could also serve to test the validity of this assumption as well as the general accuracy of our previous analytic treatment.
To make the paper self consistent, in the next section we repeat some of the results from our previous paper [18]. In particular, we describe the NC deformation by the angular twist and its consequences on a massive, charged scalar field propagating in the fixed (undeformed) RN background. Finally, we derive the NC master equation for the scalar field. This equation is our starting point in Section 3. There we investigate the QNM spectrum of the massive, charged scalar field by the WKB method. Our analysis is analytic and it is therefore limited by an approximation l + 1 qQ, where l is the orbital momentum of the scalar perturbation, Q is the black hole charge and q is the charge of the scalar perturbation. After this warm-up, in Section 4 we attack the problem of finding the QNM spectrum with the well defined numerical method, the continued fraction method. In Section 5 we discuss our results and give some final remarks. Details of a few cumbersome calculations are given in the Appendix.

NC scalar field in the RN background
A noncommutative deformation of space-time can be introduced in different ways [19]. We follow the twist approach [20] and deform the Poincaré algebra to a twisted Poincaré algebra. The algebra relations remain the same, while the comultiplication changes. This change is relevant for multiparticle states [21]. One of the advantages of the twist deformation is that it induces a deformed differential calculus in a well defined way. In particular, it introduces a deformed product of functions, the -product.
Our goal is to study the motion of a massive, charged scalar field in the geometry of the RN black hole described by the metric Here M is the mass of the RN black hole, while Q is the charge of the RN black hole. This problem requires a model of NC gravity coupled with a NC scalar and a NC electromagnetic field. A NC deformation of gravity is not easy to construct, see [22] and references therein. Therefore, we follow a semiclassical approach based on the assumption that the geometry (gravitational field) is classical (it is not deformed by noncommutativity), while the scalar field propagating in the RN background "feels" the effects of space-time noncommutativity. In order to realize this approach, we choose a Killing twist. The twist is given by The NC deformation is controlled by the small deformation parameter a ≈ l N C . Vector fields X 1 = ∂ t , X 2 = ∂ ϕ are commuting vector fields, [X 1 , X 2 ] = 0, therefore the twist (2.2) is an Abelian twist [23]. We call (2.2) an "angular twist" because the vector field X 2 = ∂ ϕ is a generator of rotations around z-axis. The vector fields X 1 and X 2 are two Killing vectors for the metric (2.1) and that is why the twist (2.2) is called a Killing twist. In particular, the twist (2.2) does not act on the RN metric and it does not act on functions of the RN metric. In this way we ensure that the geometry remains undeformed.
The twist (2.2) defines the -product of functions (fields): Using this product, one can write an action for a massive, charged scalar field in the fixed RN background as The scalar fieldφ is a complex charged scalar field with mass µ and charge q. It transforms in the fundamental representation of the noncommutative U (1) gauge transformations. Therefore its covariant derivative is defined as with NC U (1) gauge fieldÂ µ . Note that -products in √ −g g αβ g µν can all be removed since the twist (2.2) does not act on the metric tensor (2.1).
One can check that the action (2.4) is invariant under the infinitesimal U (1) gauge transformations defined in the following way: with the NC gauge parameterΛ. Our approach is perturbative and we have to expand the action (2.4) up to first order in the deformation parameter a. To do that we expand the -products in (2.4) and use the Seiberg-Witten (SW) map. SW map enables to express NC variables as functions of the corresponding commutative variables. In this way, the problem of charge quantization in U (1) gauge theory does not exist. In the case of NC Yang-Mills theories, SW map guarantees that the number of degrees of freedom in the NC theory is the same as in the corresponding commutative theory. That is, no new degrees of freedom are introduced.
Using the SW-map NC fields can be expressed as function of corresponding commutative fields and can be expanded in orders of the deformation parameter a. Expansions for an arbitrary Abelian twist deformation are known to all orders [24]. Applying these results to the twist (2.2), expansions of fields up to first order in the deformation parameter a follow. They are given by: To have a more compact equations and keep the track of the NC deformation, we introduce the 4 × 4 antisymmetric matrix θ αβ with the only non-zero components θ tϕ = −θ ϕt = a. It is also important to note that the coupling constant q between fields φ and A µ , the charge of φ, is included into A µ , namely A µ = qA µ . Using the SW-map solutions and expanding the -products in (2.4) we find the action up to first order in the deformation parameter a. It is given by The equation of motion for the field φ is obtained by varying this action with respect to φ and it is given by with Γ λ µν being the Christoffel symbols corresponding to the metric (2.1). The RN background also fixes the U (1) gauge field A µ to This is the electromagnetic potential of the point-like charge Q located at r = 0. Consequently, the only non-zero component of the field strength tensor F µν is the radial electric field F r0 = qQ r 2 .
(2.11) Furthermore, since the only non-zero components of the NC deformation parameter θ αβ are θ tϕ = −θ ϕt = a, the equation of motion (2.9) simplifies to where ∆ is the usual Laplace operator. We also introduced To simplify the notation, we will set G = 1 in the following.
In order to solve this equation we assume an ansatz [18] φ lm (t, r, θ, ϕ) = R lm (r)e −iωt Y m l (θ, ϕ) (2.14) with spherical harmonics Y m l (θ, ϕ). Inserting (2.14) into (2.12) leads to an equation for the radial function R lm (r) We note that the first line of this equation, which describes the system without deformation, corresponds to the equation for the radial function R lm analyzed in [25,26]. The NC contribution in (2.15) vanishes for a neutral scalar field, that is for q = 0, while it does not depend on the scalar field mass µ.

QNM spectrum: WKB analysis
Quasinormal modes are a particular solution of equation (2.15). They are specified by the following boundary conditions: purely incoming at the horizon and purely outgoing in the infinity. We mentioned in the Introduction that there are different ways to solve equation (2.15) and find the corresponding QNM spectrum. To warm up, we start with the WKB approach and present an analytic solution for the QNM frequencies, valid for a specific range of parameters. Later on, in Section 4 we move to a more general numerical method, the continued fraction method.

A modified tortoise coordinate
The starting point of the WKB method is a Schrödinger type equation where ψ = rR and Note that, like in all our calculations, this equation is valid up to first order in the deformation parameter a.
It is interesting to see what is the explicit form of the tortoise coordinate defined by (3.18). In order to see this, we note that up to first order in a the relation (3.18) can equivalently be written as clearly separating a required transformation into two parts. The first part is the standard Reissner-Nordström tortoise coordinate y (0) = r RN * while the second part represents the term coming exclusively from the NC deformation.
Integrating (3.21) we find where y (0) is the standard tortoise coordinate for the Reissner-Nordström metric given by We repeat once again that the result (3.22) is valid up to first order in the parameter a. The corrections in (3.22) are such that they do not change the position of the event horizons, which coincides with the analysis of the Hawking radiation in a semiclassical tunneling formalism [27,18]. The coordinate y (3.22) has the standard properties of a tortoise coordinate. As r → +∞, y → +∞, and as r → r + , y → −∞. Moreover, a brief inspection of the potential (3.20) shows that as r → +∞, the potential V tends to V → ω 2 − µ 2 . Similarly, as r → r + , the potential V approaches the value V → ω − qQ r + 2 1 − 2iam qQ r + . These two limiting values completely agree with the behavior of the effective potential found in [28] for the massive charged scalar field in the Kerr-Newman background (after putting the black hole angular momentum of the reference [28] to zero and taking the commutative limit of our result).

QNM spectrum
The WKB method is based on the similarity between the equation governing the behavior of a black hole perturbation (scalar in our case) and the Schrödinger equation in the case of a potential barrier [11,12,13]. The condition for QNMs is obtained by matching two WKB solutions on each side of the potential barrier given by −V to the solution inside the barrier, with the matching done simultaneously across both of the turning points. The matching procedure [8] leads to the following condition for the QNM frequencies: where n = 0, 1, 2, 3, .... This condition clearly involves finding an extremal value of the potential V and the value of the curvature at the extremal point [11,12]. In addition, the correction terms Λ j depend on the value of the effective potential and its derivatives (up to j-th order) in the maximum. The explicit form of the WKB corrections Λ 2 and Λ 3 can be found in [12] and of Λ 4 , Λ 5 , Λ 6 in [29]. Here, we will work in the 1st WKB order. Due to the additional terms in (3.20) induced by noncommutativity, the higher WKB corrections are very cumbersome to calculate. However, these corrections are small and they will not influence the leading order behavior. Therefore, we postpone the analysis of the higher WKB corrections for our future work. In order to proceed we follow [26] and introduce new dimensionless variables In terms of these variables the effective potential can be written as It is obvious from these definitions that H(r + ) < l(l + 1) + 1. Note that we discarded all terms quadratic in x other than the first parabolic term ∼ (x + Ω) 2 as non physical. Indeed, if we keep all quadratic terms, the resulting effective potential does not describe a realistic problem at hand. In particular, within our approximation the parabolic term ∼ (x + Ω) 2 has a peak immediately next to the horizon. Retaining other terms quadratic in x would only move the peak of the potential far away from the event horizon and this is physically undesirable. That is why, of all terms quadratic in x, we keep only the first parabolic term ∼ (x + Ω) 2 in the effective potential. Moreover, retaining terms of the order O(ax 2 ) in the potential function would only produce corrections to the position of its extremal point that are of the second order within our approximation and thus can be neglected.
As far as the effective potential of the Reissner-Nordström black hole due to charged scalar perturbations is concerned, the analysis carried out in reference [26] has been focused on the regime qQ l + 1 with a purpose of being able to deal with the condition (3.24) analytically in that particular case. For the same reason we also restrict our analysis to this regime. Putting this together with the just established relation H(r + ) < l(l + 1) + 1, it is clear that this condition translates into H(r + ) In addition, we also assume that It is reasonable to expect the QNM frequencies to be centered around the classical result ω = qQ r + . Therefore, it is natural to assume that the quantity Ω is also very small. As in the reference [26], we assume here that it is of the same order as H(r + ) q 2 Q 2 and µ 2 r 2 + q 2 Q 2 . Furthermore, our analysis includes a small NC parameter a and we assume that it is also of the same order as the former small quantities. To be more precise, we have to assume that la qQ r + 1, where l is the orbital momentum of the perturbation.
To summarize, our approximation can be written as: 1. The position of the extremum of the effective potential is determined by the condition dV dx x 0 = 0, leading to It is easily seen that within our approximation, the value of x 0 = r 0 −r + r 0 is very small, ensuing that the peak of the potential is close to the event horizon, as it should be.
The value of r 0 can be obtained from the relation x 0 = 1 − r + r 0 , leading to Up to first order in H(r + ) q 2 Q 2 , , Ω and la qQ r + this equation can be inverted to give In the same way one finds the remaining quantities that will be required in the analysis, see the relation (3.30): Once again, note that all these quantities have been written only up to first order in H q 2 Q 2 , µ 2 q 2 Q 2 , Ω and la qQ r + .
Since at the extremal point x 0 the first derivative dV dx vanishes, the QNM condition (3.24) implies Knowing that dr dy = f (r) 1 + iam qQ r and inserting the above expressions for dx dr r=r 0 , V 0 and f (r 0 ) into (3.30) leads to a coupled system of equations where Ω has been split in its real and imaginary part, Ω = Ω R + iΩ I . The solution of this system of equations is given by Using (3.25) one can express the real and the imaginary part of the QNM frequency ω. For completeness, in Figure 1 we plot the dependence of Re ω and Im ω on qQ in the case of µ = 0.05, l = 100 and a = 0.000001. It is obvious from equation (3.34) that there will be a splitting of frequencies for different values of the projection of angular momentum m. This result is in agreement with our previous findings in [18]. The frequency splittings ω ± = ω(m = ±100) − ω(m = 0) are plotted in Figure 2.
The straightforward comparison of the results obtained by the WKB method with our previous results in [18] obtained analytically in the near-extremal limit in not possible. A well known feature of the WKB method is that it gives better results for higher values of the angular momentum l. Therefore we plotted our results for the case l = 100. On the other hand, in [18] we analyzed only l = 1 and l = 2 cases. However, the qualitative comparison shows the frequency splitting due to noncommutativity in both cases and confirms our conclusions about the effects of noncommutative deformation on the QNM spectrum. Notice that, unlike in [18], the splitting in Im ω is already visible in Figure 1. Namely, we plotted three cases m = −100, m = 0 and m = 100. Since the NC corrections are proportional to am, it is obvious that the effect of noncommutativity will be larger for bigger m. The relative splitting can be estimated from the Figures 1-4 as δ Re ∼ Re ω + Re ω ∼ 10 −4 . In the similar way, for the imaginary part we estimate δ Im ∼ Im ω + Im ω ∼ 10 −2 . The splitting is obviously bigger then the splitting found in [18] and this enhancement is due to the high angular momentum l = 100. Im Ω, m 100 Im Ω, m 0 Im Ω, m 100  Finally, to compare our results with the commutative results obtained in [26], we have to make one additional approximation. Namely, we assume that 1, and at the same time being of the same order of magnitude as H(r + ) q 2 Q 2 , µ 2 r 2 + q 2 Q 2 , Ω and al qQ r + . This final approximation amounts to requiring q 2 Q 2 l 4 . Taken together, our approximations select the parameter range to be l qQ l 2 .
Solving the above system of equations for Ω R and Ω I then leads respectively to and Within this approximation, there is no contribution of noncommutativity in both Ω R and Ω I . The terms in the second line in (3.35) and (3.36) are second order "small" being a product 1 of q 2 Q 2 (H(r + )+µ 2 r 2 + ) 2 and am 1 r + in (3.35) and 1 (H(r + )+µ 2 r 2 + ) 2 and am qQ r + in (3.36). To compare the obtained results with the results in [26], we have to set the mass of the field µ = 0 and the noncommutativity parameter a = 0. After this, we indeed find the agreement with (19) and (20) in [26]. However, the results in [26] are zeroth order in small variables q 2 Q 2 H 2 , H q 2 Q 2 , Ω, while our results also contain the first order corrections in these variables.

Continued fraction method
So far we used two different methods to calculate the QNM spectrum of NC scalar perturbation of the RN geometry. Both methods have their limitations. The WKB method works well for high values of the orbital momentum l. The analytic treatment in [18] was possible only in the near-extremal limit of the RN geometry. 1 Note that qQ is always bigger then 1, since l qQ l 2 . According to our approximation al qQ r + 1, therefore al 1 r + 1 is also valid. In the same way q 2 Q 2 (H(r + )+µ 2 r 2 + ) 2 1 implies 1 (H(r + )+µ 2 r 2 + ) 2

1.
In this section we implement the continued fraction method [30,31] to determine the QNM spectrum of a massive charged scalar field around the RN black hole in the presence of the noncommutative deformation of space-time. This method is less restrictive then the previous two and we expect to find results for a wider range of parameters. We mention that the analysis of the undeformed (commutative) (un)charged scalar and Dirac QNM spectrum in the RN background by the continued fraction method can be found in [28,32,33]. To our knowledge, this is the first time that the continued fraction method is applied also to the NC deformations of the commutative case.
We start by looking for the asymptotic form of the QNM spectrum in the spatial infinity r → ∞ and near the horizon r → r + . The asymptotic form can be obtained by analytically solving equation (2.15) in the asymptotic limits of the spatial infinity r → ∞ and the event horizon r → r + and imposing the QNM boundary condition of purely incoming waves on the horizon and purely outgoing waves in the infinity.
In the r → ∞ limit equation (2.15) reduces to The solution to the equation (4.37) is given by Note that, in the limit r → ∞, the parameter Ω is fixed by the leading order solution to Ω 2 = ω 2 − µ 2 . Analogously, in the near horizon limit r → r + , the associated equation of motion (3.16) reduces to . (4.41) Here Z out and Z in are the amplitudes of the outgoing and ingoing waves, respectively, and they do not depend on r (or y). In the special case of a massless scalar field µ = 0 and a vanishing space-time deformation a = 0, these asymptotic solutions reduce to the asymptotic solutions of [32]. Equation (2.15) has an irregular singularity at r = +∞ and three regular singularities at r = 0, r = r − and r = r + . To implement Leaver's method we expand the solution in terms of powers series around r = r + . Then the radial part of the scalar field looks as (4.42) The parameters δ and are determined by demanding that the solution (4.42) satisfies the boundary conditions (4.41) at the horizon and in the infinity. From the general form (4.42) of the solution, it is clear that as r → ∞, the dominant behavior is determined by the term r e iΩr . Likewise, as r → r + , the dominant behavior of (4.42) is given by the term (r − r + ) δ . On the other side, we insert the expression (3.22) for the tortoise coordinate into (4.38) and (4.40) and then compare the resulting expressions with the formerly deduced asymptotic of (4.42). From this comparison the parameters and δ follow immediately and they are given by It is worth noting that these parameters are the same as the corresponding ones in the reference [33] and thus are not affected by the noncommutative space-time deformation.

Recurrence relations
Now we insert the power series solution (4.42) together with (4.43) into equation (2.15).
In this way we obtain the recurrence relations for the coefficients a n . The calculation leading to the recurrence relations is long and it is presented in the Appendix A. The result is the 6-term recurrence relation A n a n+1 + B n a n + C n a n−1 + D n a n−2 + E n a n−3 + F n a n−4 = 0, n 4 A 3 a 4 + B 3 a 3 + C 3 a 2 + D 3 a 1 + E 3 a 0 = 0, n = 3 The coefficients A n , B n , C n , D n , E n and F n are given by The coefficients α n , β n , γ n are Let us clarify these relations. The first relation in (4.44), for n 4, is a general 6-term recurrence relation. The remaining four relations are the indicial equations relating the lowest order coefficients a n in the general expansion (4.42). They may be thought as boundary conditions for the first relation in (4.44). The presence of NC deformation, through the terms linear in the NC parameter a, induces the 6-term recurrence relation. In the commutative case and for the non-extremal RN background the 3-term recurrence relation is obtained [32,33].
To compare the commutative limit of our result (4.44) with the results in the commutative case, we have to go back to the equation (A.8). There we impose the commutative limit a → 0 and divide by (r + −zr − ) 3 the remaining terms. The obtained equation then results in the 3-term recurrence relation of [32,33], α n a n+1 + β n a n + γ n a n−1 = 0, α 0 a 1 + β 0 a 0 = 0, (4.47) with the coefficients α n , β n and γ n given in (4.46). Note that the noncommutativity does not influence the parameters α n , β n and γ n .
Having the recurrence relations that involve more than 3 expansion coefficients a n , as we do have here, we cannot straightforwardly apply the usual method for solving the recurrence relations [34]. Instead, we should first use the Gauss elimination method to gradually reduce the initial recurrence relation from the 6-term recurrence relation to a 3-term recurrence relation. In our case the Gauss elimination method needs to be applied 3 times in a row. The details of this calculation are presented in the Appendix B. The final result is the 3-term recurrence relation n a n+1 + B (3) n a n + C (3) n a n−1 = 0, A In order to solve the 3-term recurrence relation (4.48), we start with the following observation. Since r + is a regular singular point, the general expansion (4.42) converges for r + ≤ r < ∞. Demanding convergence at r = ∞ implies that the sum n a n also converges. Therefore, defining the quantity R n by R n = − a n+1 an , we see that the infinite series (4.42) will converge if R n decreases sufficiently fast with the increase of n. Moreover, from the first relation in (4.48) we see that the coefficients a n must satisfy the recurrence relation , that is a n a n−1 (4.49) This relation, in combination with the second relation in (4.48) leads to the infinite continued fraction equation (4.50) The quantity R n may be interpreted as being the remaining part of the infinite continued fraction. It is supposed to decrease relatively fast as n grows. Hence, the convergence of the series (4.42) is ensured if the coefficients a n satisfy the equation (4.50) and R n decreases with increasing n. The solution to this infinite continued fraction equation gives the QNM frequencies. The continued fraction relation (4.50) can be inverted any number of times. Numerically, the n-th QNM frequency is defined to be the most stable root of the n-th inversion of the continued fraction relation (4.50), The fundamental mode is obtained as the most stable root of the equation (4.50).

Numerical results
Now we use the root finding algorithm (Wolfram Mathematica) to determine the QNM spectrum from equations (4.50) and (4.51). The NC parameter a is fixed to a = 0.01. We present our results graphically. The dependence of the fundamental QNM frequency ω = Re ω + iIm ω on the charge 2 of the scalar field qQ is shown in Figure 3. The remaining parameters are fixed as follows: µ = 0.05, l = 1 and M = 1 in accordance with [32,33]. The extremality of the RN black hole is controlled by the value of Q M . We present results for ten different cases, with Q M varying from 0.01 to 0.999999. One can see immediately from Figure 3 that, as the electromagnetic interaction increases (qQ increases), the imaginary part of the quasinormal frequency approaches a constant value, while the real part seems to grow linearly. In addition, as qQ decreases, the real part of the fundamental frequency Re ω approaches zero at some critical value of the electromagnetic coupling. This behavior is similar to the behavior in the commutative case [32]. As this critical value is approached, the continued fraction method seems to converge slower. The imaginary part of the fundamental frequency Im ω changes its behavior as the extremal limit is approached, at some point between Q M = 0.8 and Q M = 0.9. As the extremal limit Q M = 1 is approached, the imaginary part of the fundamental frequency becomes smaller, while the real part approaches qQ r + ≈ q.
The existence of the modes with arbitrarily small imaginary part in the nearextremal limit could be an artifact of the complicated continued fraction equation, as suggested in [32]. Namely, Leaver's original method fails to converge for extremal black holes. More precisely, when Q M → 1, the regular singularities at r = r + and r = r − merge at r = M , becoming an irregular singularity. Therefore, it is not expected that a power series expansion around r = r + will have a nonzero radius of convergence. Onozawa et al. [35] have proposed a modification of Leaver's method to deal with such type of equations and have successfully applied it to uncharged fields around an extremal RN black hole. The obtained results are in very good agreement with the ones obtained for nearly extremal black holes using Leaver's original method. In our previous work [18], we treated equation (2.15) analytically in the near extremal limit and we verified the existence of the modes with arbitrary small imaginary parts. In fact, the QNMs obtained in [18] in the near extremal limit are in excellent agreement with the ones obtained through the continued fraction method in the limit Q M → 1.  [18]. Finally, let us comment on the effects of noncommutative deformation in our model. It is clear from equations (4.44) and (4.45) that the QNM frequencies depend on the value of am, where a controls the NC deformation and m = −l, −l + 1, ..., l is the projection of the orbital momentum l. Therefore, we expect to find a frequency splitting for different values of m. To make this observation more explicit, in Figures 5 and 6 we show the frequency splitting as a function of the scalar field charge q ( Figure 5) and of the scalar field mass µ ( Figure 6). The splitting can be graphically presented for an arbitrary value of Q M , taking into account comments on the near extremal limit from the previous paragrph. In the example shown in Figures 5 and 6 we fixed the extremality to Q M = 0.5, that is we show splitting in the non-extremal case. The splitting is defined as ω ± = ω(m = ±1) − ω(m = 0).
The relative frequency splitting can be estimated from Figures 3 and 5 and Figures  4 and 6. For example, in the case of fixed mass µ = 0.05 and Q M = 0.5 we find δ Re ∼ Re ω + Re ω ∼ 10 −4 and δ Im ∼ Im ω + Im ω ∼ 10 −4 . One can check that in the case of fixed charge qQ = 1, the splitting is approximately of the same order. It is important to stress that in the limit a → 0 (and µ → 0) our results reduce to the results presented in [32,33,28], meaning that we have been able to reproduce in full detail the same graph profiles as in these references. In particular, for the uncharged massive scalar field in the absence of noncommutativity and for Q M = √ 0.1 ≈ 0.316 we checked that the curves showing the dependence of Re ω and Im ω on the mass µ are in a perfect agreement with the corresponding curves in reference [33], see the graphs on Figures 2c) and 2d) respectively, specified for the case of the scalar charge s = 0.1. In Figures 3 and 4 the dark blue lines correspond to this particular extremality Q M = 0.316. However, in this case, unlike in [33] the scalar field is charged and in addition there are effects of noncommutativity. Another example which may be interesting to confront with our results is the fundamental mode of the massless scalar field in a non-rotational Kerr-Newman background which was studied in reference [28], see Figure 3 in [28] for the case r − = 0.95r + . This case roughly corresponds to the extremality ratio of Q/M = 0.9999, which is represented on our Figure 3 by the dark green line. A close inspection of two graph profiles shows a high level of agreement. In particular, the imaginary parts in both cases saturate roughly at the same value and both have the minimum lying on the vertical axis. The small discrepancy between the two profiles may be attributed to the fact that the mass of the scalar probe in these two cases is not the same, as well as to a small effect caused by noncommutativity.

Discussion and final remarks
In this paper we studied the QNM spectrum of the NC scalar perturbation in the nonextremal RN background in more depth. In our previous work [18] we did the analytic analysis limited to the near-extremal geometry. Here we overcome this limitation by using the WKB and the continued fraction methods. The WKB method is well defined for large values of the orbital momentum l, and due to our additional approximations is limited to the parameter region l + 1 qQ. In this approximation we solved the WKB condition analytically and plotted our results in Figures 1 and 2. If one makes the additional approximation l qQ l 2 , the massless, commutative limit of our results (3.35) and (3.36) agrees with the results obtained in [26].
Looking for a less restrictive method, we then turned to the continued fraction method. To our knowledge, this is the first time the continued fraction method is applied to the noncommutative QNM spectrum problem. The NC deformation induces the 6-term recurrence relation (4.44) and we use the Gauss elimination procedure to reduce it to the 3-term recurrence relation (4.48). We used the root finding algorithm to solve the 3-term recurrence relation and obtain the fundamental QNM frequency. The results we obtain are plotted in Figures 3-6 and are in good qualitative agreement with our previous results in [18]. We discuss some of the properties of the QNM spectrum in Section 4.2. Here, for completeness, we also show the results for l = 2. In Figure 7 we plot the dependence of the fundamental frequency ω on the scalar field charge q, while in Figure 8 we plot the corresponding frequency splitting defined as ω ±± = ω(m = ±2) − ω(m = 0).  Finally, to compare the results obtained by the continued fraction method with those obtained in the WKB approximation, we present on Figure 9 the real and the imaginary part of the fundamental frequency ω in the case of l = 100 and the NC parameter fixed at a = 0.000001. The frequency splittings are presented in Figure 10. Comparing with Figures 1-2, we notice a qualitative agreement between these two methods. Quantitative results are a bit different, but that was expected. WKB is an approximative method and we made additional approximations in order to obtain an analytic result for the frequencies (3.36).
As we mentioned before, this is the first time the continued fraction method is applied to a NC QNM spectrum problem. In our future work, we plan to investigate spinor and vector NC QNM spectrum. Of course, we are interested to go beyond the semiclassical analysis used in this paper. To do that, we need a full NC gravity action. We hope that we can make some progress in this direction using the model developed in [22].
We have to keep in mind that the complete radial solution must satisfy the boundary conditions (4.41) on the horizon and at far infinity. The change of coordinates from r to z is accompanied with the transformation of derivatives Inserting further the expansion (A.3) into (A.1), and using the decomposition r = ∞ n=0 a n (n + δ)z n+δ−1 ∞ n=0 a n z n+δ + 1 (1 − z) ∞ n=0 a n (n + δ)z n+δ−1 ∞ n=0 a n z n+δ = 0.
The coefficients A n , B n , C n , D n , E n , F n that appear here are given in (4.45). Moreover, it is easy to see that the expressions multiplying a 0 , a 1 , a 2 , a 3 , a 4 that appear in front of the lowest order powers z δ , z δ+1 , z δ+2 and z δ+3 may be recognized as the particular cases of the coefficients A n , B n , C n , D n , E n , F n . Indeed, these expressions may be written in terms of A i , B i , C i , D i , E i , i = 0, 1, 2, 3. Thus, the sequence of recurrence relations following from (A.9) can be written as A n a n+1 + B n a n + C n a n−1 + D n a n−2 + E n a n−3 + F n a n−4 = 0, A 3 a 4 + B 3 a 3 + C 3 a 2 + D 3 a 1 + E 3 a 0 = 0, A 2 a 3 + B 2 a 2 + C 2 a 1 + D 2 a 0 = 0,

B Gauss elimination procedure
In the following we explain the required steps in a Gauss elimination procedure applied to our particular model. We need to reduce the 6-term recurence relation (4.44) to a 3-term recurence relation.
To begin with, we define the coefficients of the zeroth level A n , j = 1, 2, 3 as the coefficients that appear in the (6 − j)-term recurence relation, obtained after the jth Gauss elimination.
Applying the first Gauss elimination to (4.44) we find the 5-term recurrence relation A (1) n a n+1 + B (1) n a n + C (1) n a n−1 + D (1) n a n−2 + E (1) n a n−3 = 0, A 14) The application of the second Gauss elimination to the recurrence relation (B.12) yields the following 4-term recurrence equation n a n+1 + B (2) n a n + C (2) n a n−1 + D (2) n a n−2 = 0, A 1 a 2 + B n , B n , C n , that we have been searching for, are given by and for n = 0, 1,