American Journal of Computational Mathematics
Vol.2 No.2(2012), Article ID:20069,11 pages DOI:10.4236/ajcm.2012.22013

Solving Systems of Transcendental Equations Involving the Heun Functions

Plamen P. Fiziev1,2, Denitsa R. Staicova1

1Department of Theoretical Physics, Sofia University “St. Kliment Ohridski”, Sofia, Bulgaria

2B-Layer Task Force, JINR, Moscow, Russia

Email: {fiziev, dstaicova},

Received March 7, 2012; revised April 8, 2012; accepted April 16, 2012

Keywords: Root-Finding Algorithm; Müller Algorithm; Two-Dimensional Müller Algorithm; Regge-Wheeler Equation; Quasinormal Modes; Teukolsky Master Equation


The Heun functions have wide application in modern physics and are expected to succeed the hypergeometrical functions in the physical problems of the 21st century. The numerical work with those functions, however, is complicated and requires filling the gaps in the theory of the Heun functions and also, creating new algorithms able to work with them efficiently. We propose a new algorithm for solving a system of two nonlinear transcendental equations with two complex variables based on the Müller algorithm. The new algorithm is particularly useful in systems featuring the Heun functions and for them, the new algorithm gives distinctly better results than Newton’s and Broyden’s methods. As an example for its application in physics, the new algorithm was used to find the quasi-normal modes (QNM) of Schwarzschild black hole described by the Regge-Wheeler equation. The numerical results obtained by our method are compared with the already published QNM frequencies and are found to coincide to a great extent with them. Also discussed are the QNM of the Kerr black hole, described by the Teukolsky Master equation.

1. Introduction

The Heun functions appear with increasing frequency in modern physics. For example, they arise in the Schrö- dinger equation with anharmonic potential, in water molecule, in the Stark effect, in different quantum phenomena related with repulsion and attraction of levels, in the theory of lunar motion, in gravitational physics of scalar, spinor, electromagnetic and gravitational waves in Schwarzschild and Kerr metric, in crystalline materials, in three-dimensional waves in atmosphere, in Bethe anzatz systems, in Collogero-Moser-Sutherland systems, e.t.c., just to mention a few. Because of the wide range of their applications ([1,2])—from quantum mechanics to astrophysics, from lattice systems to economics—they can be considered as the 21st century successors of the hypergeometric functions encountered in some simple physical problems of 20th century.

It is not hard to explain this situation. In natural sciences, in particular, in physics we usually study the different phenomena starting from some equilibrium state. Then we study small deviations from it in linear approximation, and at the end, going far away from the equilibrium we are forced to take into account nonlinear phenomena. It is well known that to describe the wave processes (like those in quantum mechanics), related with some linear phenomenon in classical physics (like classical mechanics), we have to use hypergeometric functions. Therefore these functions were well studied in 19th and 20th centuries and today one can find the corresponding codes in all good computer packages. According to the theorem by Slavyanov [2], if we study nonlinear classical phenomena, described by elliptic functions, or even by the solutions of any of Painlevé type equations, the corresponding wave problems can be solved exactly in terms of the Heun functions. Since the Painlevé equations can be considered as Hamilton ones for a very large class of nonlinear classical problems, one can expect a fast increase in the applications of the Heun functions in physics and other natural sciences of 21st century.

Their mathematical complexity, however, makes working with them a significant challenge both analyticcally and numerically. The Heun functions are unique local Frobenius solutions of a second-order linear ordinary differential equation of the Fuchsian type [2-5] which in the general case have 4 regular singular points. Two or more of those regular singularities can coalesce into an irregular singularity leading differential equations of the confluent type and their solutions: confluent Heun function, biconfluent Heun function, double confluent Heun function and triconfluent Heun function. The Heun functions generalize the hypergeometric function (and also include the Lamé function, Mathieu function and the spheroidal wave functions [2,5]) and some of their uses can be found in [2] and also in the more recent [6]. Clearly, the Heun functions will be encountered more and more in modern physics, hence, there is a need of better understanding of those functions and new, more adequate algorithms working with them.

Despite the growing number of articles which use equations from the Heun type and their solutions, the theory of those functions is far from complete. There are some analytical works on the Heun functions, but they were largely neglected until recently. Some recent progress can be found in [7], but as a whole there are many gaps in our knowledge of those functions. Particularly, the connection problem for the Heun functions is not solved—one cannot connect two local solutions at different singular points using known constant coefficients [2]. Another example of a serious gap in the general theory of the Heun functions in general is the absence of integral representations analogous to the one for hypergeometric functions. According to Whittaker’s hypothesis, the Heun’s functions are the simplest class of special functions for which no representations in form of contour integrals of elementary functions exists. In some particular cases are known integral representations in form of double integrals of elementary functions, but the general case is an open problem.

Numerically, the only software package currently able to work with the Heun functions is Maple. Alternative ways for evaluations of those functions do not exist (to the best of our knowledge) and there are no known projects aiming to change this situation, an admittedly immense task by itself. This means that the use of the Heun functions is limited to the routines hidden in the kernel of Maple, which the user cannot change or improve—a situation that makes understanding the numerical problems or avoiding them adequately very difficult. On the positive side, those routines were found by the team to work well enough in many cases (see, for example, the match between theory and numerical results in [8], as well as the other applications of those functions in [9,10]). Yet, there are some peculiarities—there are values of the parameters where the routines break down leading to infinities or to numerical errors. The situation with the derivatives of the Heun functions in Maple is even worst—for some values they simply do not work, for example outside the domain, where their precision is much lower than that of the Heun function itself. Also, in some cases there are no convenient powerseries representations and then the Heun functions are evaluated in Maple using numerical integration. Therefore the procedure goes slowly in the complex domain (compared to the hypergeometric function) which means that the convergence of the root-finding algorithm is essential when one solves equations including Heun’s functions.

Despite all the numerical set-backs, the Heun functions offer many opportunities to modern physics. They occur in the problem of quasi-normal modes (QNM) of rotating and nonrotating black holes, which is to some extent the gravity analogue of the problem of the hydrogen atom. Finding the QNMs is critical to understanding observational data from gravitational wave detectors and proving or refuting the black holes existence ([11,12] and also [8-10]). In this case, one has to solve a two-dimensional connected spectral problem with two complex equations in each of which one encounters the confluent Heun functions. The analytical theory of the confluent Heun function is more developed than that of the other types of Heun functions, but still many unknowns remain. Again, the evaluation of the derivative of the confluent Heun function is problematic in Maple, which makes Newton’s root-finding algorithm ([13,14]) unusable. Broyden’s algorithm ([15]) generally works, but it is slowly convergent even close to a root (see [16]). It is clear that we need a novel algorithm, that will offer quicker convergence than Broyden’s algorithm, but without relying on derivatives.

To solve this problem, in the case of a system of two equations in two variables, our team developed a twodimensional generalization of the Müller algorithm. The one-dimensional Müller algorithm ([17]) is a quadratic generalization of the secant method, that works well in the case of a complex function of one variable. It has very good convergence for a large class of functions (~1.84) and it is very efficient when the starting point (the initial guess) is close to a root (for applications see [18] and [8]). It is also well convergent when working with special transcendental functions such as the Heun functions. Most importantly, this algorithm does not use derivatives, which is key for our work.

The two-dimensional Müller algorithm seems to inherit some of the advantages of its one-dimensional counterpart like good convergence and usability on large class of functions as our tests show. The new algorithm was used to solve the QNM problem in the case of a Schwarzschild black hole and it proved to work without significant deviations from the results published by Andersson ([19]) and Fiziev ([18]). Also, preliminary results for the QNM of the Kerr black hole are discussed and for them we also obtain a very good coincidence with published results [20].

The article is organized as follows: Section 2 reviews the one-dimensional Müller algorithm and it introduces in detail its two-dimensional generalization, in Section 3 we give two examples of the application of the new method in systems featuring the confluent Heun functions—the QNMs of rotating and nonrotating BH, and the results in both cases are analyzed in terms of precision and convergence. In Section 4 we summarize our results.

2. The Müller Algorithm

2.1. One-Dimensional Müller’s Algorithm

The one-dimensional Müller algorithm ([13,17]) is iterative method which at each step evaluates the function at three points, builds the parabola crossing those points and finds the two points where that parabola crosses the x-axis. The next iteration is the the point farthest from the initial point.

Explicitly, the one-dimensional Müller algorithm can be defined as the map, which obtains the final point in iterations by calculating for every three points, , (with the corresponding values of the function, ,), the next iteration as:


The iterations continue until where is the number of digits of precision we require. This exit-condition works independently of the actual numerical zero in use, which may vary for the confluent Heun function and thus it is the most appropriate for our numerical work.

2.2. Two-Dimensional Müller’s Algorithm

The two-dimensional Müller method comes as a natural extension of the one-dimensional Müller method.

For two complex-valued functions and we want to find such pairs of complex numbers which are solutions of the system:


where numbers the solution in use. From now on, we will omit the index, considering that we work with one arbitrary particular solution. Finding all the solutions of a system is beyond the scope of this article.

Consider the functions as twodimensional complex surfaces and

in a three-dimensional space of the complex variables 1. Normally, to solve the system, one expresses the relation from one of the equations, then by substituting it in the other equation, one solves it for and from one finds. In the general case, however, this is not possible. The idea of our code is to approximately follow that procedure by finding an approximate linear relation between the two variables and then using it to find the root of function of one variable trough the one-dimensional Müller algorithm.

To find the linear relation, at each iteration we form the plane passing trough three points of one of the functions and then the equation of the line of intersection between that plane and the plane is used as the approximate relation. This basically means that the so found is an approximate solution of one of the equations which ideally should be near the real solution in the plane. Substituting this relation in the other function, we run the one-dimensional Müller algorithm on it to fix the value of one of the variables, say. Using the value of in the first function, we again run the one-dimensional Müller algorithm on it to fix the value of the other variable—y. Alternatively one can substitute the value of directly in to obtain. This ends one iteration of the algorithm. The process repeats until one of the exit-conditions discussed below has been reached. The block-scheme of the algorithm can be seen on Figure 1.

Explicitly, the code starts by evaluating the two functions in three starting pairs of points () that ideally should be near one of the roots of the system. In our case, those three initial pairs are obtained from one starting pair to which we add and subtract certain small complex number. This artificial choice is done only in the first iteration (), afterwards we use the output of the last three iterations to form, and the respective. Thus on every iteration after the actual complex functions are evaluated only once outside of the one-dimensional Müler subroutines.

Next we construct the plane passing trough those three points for one of the functions, say by solving the linear system:

Figure 1. A block scheme of the two-dimensional Müller algorithm. λ(C1, C2, C3) is the plane with equation z = C1x + C2y + C3 that crosses trough the pairs of points (xi,yi) and the function F2 evaluated in them. The plane v is defined by the equation z = 0. The one-dimensional Müller algorithm, μ(tin, F(t)) → tfin, is applied on the function of one variable F(t) with starting point tin and final point tfin.

From it one obtains the coefficients of the plane.

This plane is intersected with the plane and the equation of the line between those two planes is the approximate relation (i.e.) of the two variables.

We substitute that relation in the first function and we start the one-dimensional Müller on that “linearized” function of only one variable, x. After some pre-determined maximal number of iterations, the exiting point is chosen for () .

Then, there are two possibilities.

Algorithm M1: one could use directly the relation to find. OrAlgorithm M2: One can substitute in the other function in order to find using again the one-dimensional Müller algorithm.

Our numerical experiments showed that both approaches lead to convergent procedure.

After are fixed, the two functions are evaluated and if the new points are not roots, the iterations continue.

The exit-strategy in the two-dimensional Müller algorithm is as follows:

1) To avoid hanging of the algorithm or its deviation from the actual root of the system, one fixes the maximal number of iterations for the one-dimensional Müller subroutine, P. From our experience, P = 4 - 10 gives best convergence.

2) The precision-condition remains in force for the one-dimensional Müller algorithm. Usually the subroutine exits, because of during the first few iterations of the two-dimensional Müller algorithm, until it gets closer to the roots and exits with.

3) The primary exit-condition of the two-dimensional Müller is fulfilled when for two consecutive pairs:. In this case, if the values of functions at those points are sufficiently small, the algorithm exits with a root.

4) To keep the two-dimensional Müller algorithm from hanging, a maximal number of iterations should be set. For the algorithm exits without fixing a root.

5) A common problem occurs when one of the functions becomes zero before the other function. A possible way out is to substitute one of the so-fixed variables, say, in the non-zero function and to to run the onedimensional Müller subroutine using the other variable—. The algorithm then exits with a possible root:.

The procedure can be fine-tuned trough change in the starting pair of points, the initial deviation or by switching the places of the functions, or even by replacing the functions with their independent linear combinations.

The numerical experiments of the application of this algorithm on systems featuring different classes of functions can be found in [16]. The tests showed that the algorithm inherits some of the advantages of the one-dimensional Müller algorithm, like the quick convergence in proximity of the root and the vast class of functions that it can work with. The major disadvantage comes from the complicated behavior of the two-dimensional complex surfaces defined by the functions which requires one to find the best combination of starting points and number of iterations in the one-dimensional Müller subroutine so that the algorithm converges to the required root (if it is known or suspected). Generally, it is hard to tell when one point is “close” to a root. In some cases, even if certain starting pair of points is close to a root in terms of some norm, using it as a starting point in the algorithm may still lead to convergence to another root or simply to require more iterations to reach the desired root than if other pair of starting points were used.

It is important to note that unlike Broyden’s algorithm and Newton’s algorithm which do not dependent on the order of the equations in the system, our two-dimensional Müller algorithm depends on the order of the equations. The numerical experiments show that while for some systems, changing the places of the equations has little or no effect on the convergence, in other cases, it slows down or completely breaks down the convergence. While such inherent asymmetry certainly is a weakness of the algorithm, there are ways around it. For example, one may alternate the places of the equations at each iteration or use their independent linear combinations . Those approaches make the algorithm more robust, but since they may cost speed, we prefer to set the order of the equations manually.

A technical disadvantage is that the whole procedure is more CPU-expensive than Newton’s method and Broyden’s method, since it generally makes more evaluations of the functions—each one-dimensional Müller makes at least 1 iteration on every step of the two-dimensional Müller, thus it makes at least 4 evaluations of each function. This is because on each iteration of the two-dimensional Müller algorithm the functions in use change and thus one cannot use previous evaluations to reduce time. Still, in some cases, as demonstrated in [16] and also below, the so-constructed algorithm is quicker or comparable to Newton’s or Broyden’s method.

3. Some Applications of the Method—QNMs of Nonrotating and Rotating Black Holes

We will work only with the confluent Heun functions, which are much better studied than the other types of Heun functions, due to their numerous physical applications. Besides their numerical implementation was used successfully in previous works by the authors. For details on the numerical testing, see [16].

The quasi-normal modes (QNMs) of a black hole (BH) are the complex frequencies that govern the latetime evolution of the perturbations of the BH metric ([11,12,21-23]), which have been extensively studied.

3.1. First Example: Non-Rotating Black Hole

First, we consider the problem of the gravitational QNMs of a nonrotating black hole so that the precision of the new method can be tested on a very well studied physical problem. The physical results in this case were published in [9], so here we will focus on the numerical details instead.

To find the QNMs, one uses the exact solutions of the Regge-Wheeler equations, describing the linearized perturbations of Schwarzschild metric, in terms of confluent Heun functions([18]). From [18], when the mass of the BH is set to, one obtains the following system of the type (1):


where is a complex frequency, is the angular momentum of the perturbation, is the angle which we set to and. HeunC is the confluent Heun function ([3]) in Maple notations. The parameter is a small variation in the phase condition (see [18]).

Using Equation (2), we run the two-dimensional Müller algorithm to find the unknown and with precision of the algorithms set to 15 digits.

From the theory, it is known that is an integer and. Comparing with the results obtained by the two-dimensional Müller algorithm, for the first root, one has, with the first different from 9 digit being the 17th. This shows that the new algorithm is capable of solving systems with one purely integer root in the pair with the expected precision.

A comparison of the new algorithm with the wellknown Newton’s and Broyden’s methods, can be found in Table 1.

Because the phase condition includes the complex argument in non-analytical way, which cannot be differentiated, this problem cannot be solved directly using Newton’s method. Broyden’s method works but with serious limitation of its precision. This happens, because one of the roots in the pair is a real integer, while the other is complex and the algorithm fixes the integer root very quickly, thus the finite differences in the Jacobian become infinity. Because of this, the algorithm is able to fix only the first 10 - 11 digits, while the other algorithms fix 14 - 15 digits. Therefore, although Broyden’s algorithm gives better times (see Table 1) than the two-dimensional Müller algorithms, its precision is much lower and for modes with big imaginary part, it cannot be increased even by raising the software floating-point number to very high

Table 1. The times needed for Broyden’s method (tB) and the two-dimensional Müller methods (tM1 and tM2) to fix a root. Note that while the precision of the former is 10 digits, the precision of the other two is 14 - 15 digits. To obtain those times, we solve the system: [F1 + F2, F1 – F2] with starting points: ω[n] + 0.01 + 0.01i, 2.1 + 0.01i, where n = 0.10.

values. Furthermore, from Table 1, one can see that the time needed for the each algorithm to exit with a root dramatically increases with. This emphasize on the importance of the convergence of the algorithm, which may become critical in physical problems where multiple roots must be found (see the second example).

The numerical results for the QNMs are summed in Table 2. In it, the QNM frequencies obtained from Sys. (2) are compared to those found by Andersson ([19]) with the phase amplitude method. Recently, those results were confirmed by Fiziev (see [18]) with the one-dimensional Müller method applied on the exact solutions of the radial equation in terms of the confluent Heun function for. To check the accurateness of the new method, we evaluate.

From the table, it is clear that in most cases, the modes obtained with the two-dimensional Müller algorithm coincide with those obtained by Andersson with more than 8 digits of precision in most cases and for modes there are 9 coinciding digits (Andersson published 9 digits of his frequencies). These results also confirm the roots for published in [18].

The mode with biggest deviation from the expected value is in Table 2 and it was already discussed in [9] (and references therein). In brief its properties are due to the branch cuts in the radial function, which also lead to non-trivial dependence of the frequencies on (where and): for ); for

and for,


Because of this, the value for in the table 2 was obtained for positive (), unlike the other modes with, which were obtained for.

The frequencies presented here are stable with precision of 6 digits at the worst and usually around 9 digits with respect to a change of in the corresponding intervals.

To confirm that the so-observed dependence of on the parameter is not a problem of the algorithm, but it is a feature of the numerical realization of the confluent Heun function, we make complex 3d plots of the radial function in use for several values of. The effect of this parameter on the branch cuts in the radial function was already discussed in [9] and [10], so on Figures 2-4 in the Appendix we only illustrate the movement of the branch cuts under the change of.

The appearance of those branch cuts represents one more complication in the work with the confluent Heun functions, since not all of them are documented.Using the parameter, however, those branch cuts can be controlled and one can obtain valuable physical results.

3.2. Second Example: Rotating Black Holes

A significantly more complicated system to solve can be found in the case of QNMs from rotating black holes. In this case, one uses the exact solutions of the Teukolsky radial and angular equations, describing the linearized electromagnetic perturbations of the Kerr metric, in terms of confluent Heun functions, as stated for the first time in full detail in [18]. The two-dimensional Müller algorithm was applied successfully in this case too and the complete results and their discussion can be found in [10]. Here, we present some details on the numerical procedures used in this case.

From [24], for the values of the parameters: s = –1, , , m = 0, a = 0.01, , one obtains:

Table 2. A list of the frequencies we obtained for the QNMs of Schwarzschild black hole compared with the numbers found by Andersson.. The first 5 frequencies (n = 0 - 4, marked with *) were obtained also by Fiziev using the confluent Heun functions and coincide with the presented here. The 8th mode, marked with **, was obtained by Leaver [25].

where is the derivative of the confluent Heun function ([3]) as defined in Maple.

For brevity, here the radial equation was rounded to only 4 digits of significance. In our numerical experiments, we used the complete system with software floating-point number set to 64, where the derivatives of the confluent Heun functions were replaced with the associate confluent Heun function according to equation (3.7) of [7]. This was done to avoid the numerical evaluation so that the peculiarities of the numerical implementation of the confluent Heun function (i.e. the use of Maple fdiff procedure) are minimized. The difference in the times needed to fix a root when is used and when it is not used is small for the modes (i.e. x) with small imaginary part (~ 15 s), but it increases with the mode number, until it becomes significant for modes with big imaginary part (for the 10th mode—in Table 3 the difference is already). This slowdown is due to the timeconsuming numerical integration and differentiation in the complex domain, needed for the evaluation of.

For that system, three pairs of starting points were used: (0.49 + 0.18i, 2.001 + 0.1i), (0.17 + 0.97i, 2.001 + 0.1i), (0.069 + 5.146i, 2.001 + 0.051i) The results can be found in Table 3. One sees that the two modifications of the two-dimensional Müller algorithm M1 and M2 are much quicker than the Broyden algorithm (tM1 ~ tM2 < tB) Newton ’s method once again cannot be used.

The supremacy of the Müller algorithms is clear and it is not isolated—it is observed for other modes or values of the parameters (for example, for). To check the precision of the method, the first two modes were compared with the already published results of electromagnetic QNMs of a Kerr black hole (see [20]) and were found to coincide with at least 9 digits of significance with them. We could not find a published value for the third mode.

This example shows that in systems featuring the confluent Heun functions and their derivatives, the two-dimensional Müller method is much better suited than the already known algorithms.

4. Conclusions

The confluent Heun function appear in many physical problems. Because of the peculiarities of those functions, however, the standard numerical root-finding algorithms do not work efficiently enough on them. Here, we presented the general idea of a method for solving a system of two complex-valued nonlinear transcendental equations with complex roots based on the one-dimensional Müller method. This method avoids some of the problems accompanying the work with the Heun functions and it is aimed to provide adequate way to deal with systems featuring those functions.

In the current article, the numerical results from the application of the new algorithm to systems from the QNM

Table 3.QNMs of Kerr BH for s = –1. R numbers the root, t and N label the time and the iterations needed for the algorithms to exit. * denotes the roots dependent on the order of the equations in M1 and M2.

physics are presented. They showed that in those cases, the new method indeed works better than the standard methods. Therefore, the new method can be readily applied to find the roots of the Regge-Wheeler equation [18], the Zerilli equation [28], the Teukolsky radial and angular equations [24], all of which are solved analytically in terms of confluent Heun functions. Using this algorithm, we were able to solve directly the problem of quasi-normal modes of a Schwarzschild ([9]) and Kerr black hole ([10]) with higher precision than that of the Broyden method. The so found solutions agree to great extent with previous published numerical results thus confirming the usefulness of the method.

The complete mathematical investigation of the proposed new method, and especially its theoretical order of convergence under proper conditions on the class of functions is still an open problem.

For other applications of the method see the recent references [29,30].

5. Acknowledgements

The authors are grateful to prof. Hans Petter Langtangen for the critical reading of the early version of the manuscript and for the useful advices.

This article was supported by the Foundation “Theoretical and Computational Physics and Astrophysics”, by the Bulgarian National Scientific Fund under contracts DO-1-872, DO-1-895, DO-02-136, and Sofia University Scientific Fund, contract 185/26.04.2010.

6. Author Contributions

P.F. gave the idea and the outline of the two-dimensional generalization of the Müller algorithm, chose the physical problem on which to test the algorithm and he supervised the project.

D.S. is responsible for the realization of the algorithm in Maple code, the testing and the optimization of the code and for the numerical results and plots presented here.

Both authors discussed the results of the tests of the algorithm, commented them and were involved in trouble-shooting of the code at all stages. The manuscript was prepared by D.S. and edited by P.F.


  1. F. W. J. Olver, D. W. Lozier, R. F. Boisvert and C. W. Clark, “NIST Handbook of Mathematical Functions,” Cambridge University Press, Cambridge, 2010.
  2. S. Y. Slavyanov and W. Lay, “Special Functions, A Unified Theory Based on Singularities,” Oxford Mathematical Monographs, Oxford, 2000.
  3. K. Heun, “Zur Theorie der Riemann'schen Functionen zweiter Ordnung mit vier Verzweigungspunkten,” Mathematische Annalen, Vol. 33, No. 2, 1889, pp. 161-179. doi:10.1007/BF01443849
  4. A. Decarreau, M. C. Dumont-Lepage, P. Maroni, A. Robert and A. Roneaux, “Formes canoniques des equations confluentes de l'equation de Heun,” Annales de la Societe Scientifique de Bruxelles, Vol. 92, 1978, p. 53.
  5. A. Decarreau, P. Maroni and A. Robert, Heun’s Differential Equations, A. Roneaux, ed., Oxford University Press, Oxford, 1995, p. 354.
  6. M. Hortacsu, “Heun Functions and Their Uses in Physics,” 2011.
  7. P. P. Fiziev, “Novel relations and New Properties of confluent Heun’s functions and Their Derivatives of arbitrary order,” Journal of Physics A: Mathematical and Theoretical, Vol. 43, No. 3, 2010, Article ID: 035203. doi:10.1088/1751-8113/43/3/035203
  8. D. Staicova and P. Fiziev, “The Spectrum of Electromagnetic Jets from Kerr Black Holes and Naked Singularities in the Teukolsky Perturbation Theory,” Astrophysics and Space Science, Vol. 332, No. 2, 2011, pp. 385-401. doi:10.1007/s10509-010-0520-x
  9. P. Fiziev and D. Staicova, “Application of the confluent Heun functions for finding the QNMs of nonrotating Black Hole,” Physical Review D, Vol. 84, No. 12, 2011, Article ID: 127502. doi:10.1103/PhysRevD.84.127502
  10. D. Staicova and P. Fiziev, “New results for Electromagnetic Quasinormal Modes of Black Holes,” 2011 arXiv:1112.0310v2
  11. S. Chandrasekhar and S. L. Detweiler, “The Quasi-Normal Modes of the Schwarzschild black hole,” Proceedings of the Royal Society A, Vol. 344, No. 1639, 1975, pp. 441-452. doi:10.1098/rspa.1975.0112
  12. S. Detweiler, “Black holes and gravitational waves. IIIThe resonant frequencies of Rotating Holes,” Astrophysical Journal, Vol. 239, 1980, pp. 292-295. doi:10.1086/158109
  13. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, “Numerical Recipes,” Cambridge University Press, Cambridge, 1992.
  14. G. E. Forsythe, M. A. Malcolm and C. B. Moler, “Computer Methods for Mathematical Computations,” Prentice Hall, Upper Saddle River, 1977.
  15. C. G. Broyden, “A Class of Methods for Solving Nonlinear Simultaneous Equations Math,” Mathematics of Computation, Vol. 19, 1965, pp. 577-593.
  16. P. Fiziev and D. Staicova, “Two-Dimensional Generalization of the Muller Root-Finding Algorithm and Its Applacations,” 2011.
  17. D. E. Müller, “A Method for Solving Algebraic Equations Using an Automatic Computer,” Mathematical Tables and Other Aids to Computation, Vol. 10, No. 5, 1956, pp. 208-215.
  18. P. P. Fiziev, “Exact Solutions of Regge-Wheeler Equation and Quasi-Normal Modes of Compact Objects,” Classical and Quantum Gravity, Vol. 23, No. 7, 2006, pp. 2447-2468. doi:10.1088/0264-9381/23/7/015
  19. N. Andersson, “A Numerically Accurate Investigation of Black-Hole Normal Modes,” Proceedings of Mathematics and Physical Sciences, Vol. 439, No. 1905, 1905, pp. 47- 58.
  20. E. Berti, V. Cardoso and C. M. Will, “On gravitationalWave Spectroscopy of massive black holes with the Space Interferometer LISA,” Physical Review D, Vol. 73, No. 6, 2006, Article ID: 064030. doi:10.1103/PhysRevD.73.064030
  21. S. Chandrasekhar, “The mathematical theory of Black Holes,” Oxford University Press, Oxford, 1983.
  22. V. Ferrari and L. Gualtieri, “Quasi-Normal Modes and Gravitational Wave Astronomy,” General Relativity and Gravitation, Vol. 40, No. 5, 2008, pp. 945-970. doi:10.1007/s10714-007-0585-1
  23. R. A. Konoplya and A. Zhidenko, “Quasinormal modes of Black Holes: From Astrophysics to String Theory,” Reviews of Modern Physics, Vol. 83, No. 3, 2011, pp. 793-836. doi:10.1103/RevModPhys.83.793
  24. P. P. Fiziev, “Classes of exact solutions to the Teukolsky master equation,” Classical and Quantum Gravity, Vol. 27, No. 13, 2010, Article ID: 135001. doi:10.1088/0264-9381/27/13/135001
  25. E. W. Leaver, “An Analytic Representation for the QuasiNormal Modes of Kerr Black Holes,” Proceedings of the Royal Society A, Vol. 402, No. 1823, 1985, pp. 285-298. doi:10.1098/rspa.1985.0119
  26. E. Berti, V. Cardoso and A. O. Starinets, “Quasinormal modes of Black Holes and Black Branes,” Classical and Quantum Gravity, Vol. 26, No. 16, 2000, Article ID: 163001. doi:10.1088/0264-9381/26/16/163001
  27. E. Berti, “Black Hole Quasinormal Modes: Hints of Quantum Gravity?” 2004.
  28. P. P. Fiziev, “Teukolsky-Starobinsky identities: A novel derivation and generalizations,” 2009.
  29. M. Cadoni and P. Pani, “Holography of Charged Dilatonic Black Branes at Finite Temperature,” Journal of High Energy Physics, 2011.
  30. P. Pani, “Applications of Perturbation Theory in Black Hole Physics,” Ph.D. Thesis, Universita’ degli Studi di Cagliari, Cagliari, 2011.


The Epsilon-Method


Figure 2. 3d plots of the function F2, the solution of the RWE, in the complex interval ω = 0.32 + 1.4i. 0.5 + 2.4i for ε = 0.05, 0, –0.05, –0.1 (the colors encode the phase of the complex function F2). The wall characteristic of the branching of the multivalued function is moved by ε either to the left (ε < 0) or to the right (ε > 0). Note that on Figures 2-4,. (a) [ε = –0.1], (b) [ ε = –0.05], (c) [ε = 0], (d) [ ε = 0.05].


Figure 3. 3d plots of the function F2 forε = –0.05, ε = 0 and ε = 0.05 in the same interval for the complex ω as in 2 scaled near the expected roots. The different ε lead to different profiles of the roots.


Figure 4. Plot of the level curves of the function for ε = 0.05, –0.05, –0.08 for ω = 0.32 + 1.4i. 0.55 + 2.4i. The movement of the branch cut due to the change of ε can be clearly seen. Note that on these plots, F2 is the radial solution of the Teukolsky Radial Equation [9], but as seen from 2, the branch cuts for this choice of parameters coincide with those of the solution of the RWE. (a) [ε = –0.08]; (b) [ ε = –0.05]; (c) [ε = 0]; (d) [ε = 0.05].


1Equivalently, we can consider four real surfaces in five-dimensional  real hyperspace, which are defined by four real functions of four real variable.