Quasinormal modes of relativistic stars and interacting fields
Abstract
The quasinormal modes of relativistic compact objects encode important information about the gravitational response associated with astrophysical phenomena. Detecting such oscillations would provide us with a unique understanding of the properties of compact stars, and may give definitive evidence for the existence of black holes. However, computing quasinormal modes in realistic astrophysical environments is challenging due to the complexity of the spacetime background and of the dynamics of the perturbations. We discuss two complementary methods for computing the quasinormal modes of spherically symmetric astrophysical systems, namely, the direct integration method and the continuedfraction method. We extend these techniques to dealing with generic coupled systems of linear equations, with the only assumption being the interaction between different fields is effectively localized within a finite region. In particular, we adapt the continuedfraction method to include cases where a series solution can be obtained only outside an effective region. As an application, we compute the polar quasinormal modes of boson stars by using the continuedfraction method for the first time. The methods discussed here can be applied to other situations in which the perturbations effectively couple only within a finite region of space.
pacs:
04.30.Db, 04.25.Nx, 04.80.Nn, 95.35.+dI Introduction
The study of the natural oscillations of physical systems is of much interest in astrophysics Berti et al. (2009); Kokkotas and Schmidt (1999). For relativistic compact objects such as stars and black holes, these natural oscillations are referred to as quasinormal modes (QNMs). QNMs depend only on the properties of the object. The radiation emitted in the collapse of a star and in the coalescence of binaries has a close relation to the QNMs of the final compact object. The composition of matter forming astrophysical objects is very important for the dynamics of the spacetime oscillations. Indeed, many matter oscillations in Newtonian theory are also necessarily present in the relativistic case Kokkotas and Schmidt (1999). Fluid modes can also be excited, for instance, by a particle moving around the compact object Kojima (1987); Pons et al. (2002), generating a large dephasing in the gravitationalwave signal Macedo et al. (2013a).
It is also of great interest to study how matter and fields around black holes behave. The presence of matter around black holes can lead to resonances, which were analyzed in detail in Ref. Barausse et al. (2014). Scalar fields around rotating Kerr black holes can generate hairy configurations Herdeiro and Radu (2014) (see also Refs. Herdeiro and Radu (2015); Brito et al. (2015a) for a recent review of black holes with scalar hair). Moreover, massive scalar fields can form quasibound states and develop scalar “clouds” around black holes Hod (2012); Benone et al. (2014); Barranco et al. (2012); Zilhão et al. (2015); Brito et al. (2015a). Notwithstanding, relativistic stars can also present nontrivial gravitating field configurations. For instance, very strong magnetic fields in neutron stars—configurations known as magnetars Duncan and Thompson (1992)—have an influence on the shape of the star Konno et al. (1999); Yazadjiev (2012). Neutron stars can also accrete dark matter, acquiring a core formed by an additional gravitating component Bertone and Fairbairn (2008); Brito et al. (2015b, 2016a). Moreover, stars formed by fundamental fields, like boson stars (BSs) Kaup (1968); Schunck and Mielke (2003); Liebling and Palenzuela (2012) and Proca stars Brito et al. (2016b), are also interesting examples of the outcomes of selfinteracting, gravitating fields.
All of the above examples share a common feature: they are formed by components (matter and fields) coupled to each other, with such a coupling being crucial for the structure of the compact object. These interactions usually have a characteristic length scale so that, at large distances, all fields decouple or only the gravitational field survives (through a powerlaw falloff).
The coupling among different fields also impacts the analysis of linear perturbation in such spacetimes. Additional components describing the perturbations of relativistic objects can enrich their oscillation spectra, generating distinctive signatures. For instance, the scalar field modes in BSs can be excited by an orbiting particle Macedo et al. (2013a). Notwithstanding, the scattering of massive scalar wave packets can carry information about the clouds around black holes Degollado and Herdeiro (2014). However, it is difficult to analyze the QNM spectrum of such systems, basically due to the fact that the coupling between the components is highly nontrivial in most cases.
Here, we present two distinct methods for computing the QNMs of relativistic systems—namely, the direct integration method and the continued fraction method—to deal with situations in which the linearized perturbations interact only in a finite region of space. The direct integration method discussed here is similar to the one presented in Ref. Pani (2013). On the other hand, we extend the continuedfraction method adopted, e.g., in Refs. Benhar et al. (1999); Nollert (1999), to deal with arbitrary coupled systems of linear equations. One of the main advantages of our improved method is that it does not need to assume a continuedfraction representation that describes the solutions in the entire domain. Instead, we only require that a continuedfraction representation exists outside an effective interaction region. This is an advantage because, in many physical situations, the metric is analytically known only outside an effective radius, so that a continuedfraction representation is available only in a portion of the full domain. We shall assume that the angular part of the perturbations can be separated such that we end up with a system of coupled second order differential equations. This separation naturally occurs in spherically symmetric configurations, and also for spinning, axisymmetric objects in the slowrotation approximation Pani (2013). Additionally, we assume that the background is stationary, so that the perturbations can be Fourier decomposed with time dependence , where is the Fourier frequency. The above assumptions are valid in a wide variety of astrophysical scenarios Berti et al. (2009); Pani (2013); Barausse et al. (2014).
The decoupling of the perturbations outside an effective radius occurs naturally in BSs. In fact, the coupling between the gravitational and scalar field perturbations decrease exponentially with the radial coordinate and become negligible beyond a certain radius, where the scalar field which composes the BSs becomes small. Therefore, as an application of our improved methods, we compute the QNMs of miniBSs through the continuedfraction technique, extending previous results obtained with less robust methods Macedo et al. (2013b). We compute the polar QNMs of miniBSs as functions of the central field density and show that monopolar modes become unstable beyond a critical density which corresponds to the maximum mass of the BS.
The remainder of this paper is organized as follows: in Sec. II we describe the direct integration and continuedfraction methods for computing the QNMs modes of compact objects surrounded by interacting fields. In Sec. III we apply the improved direct integration and continuedfraction methods to compute the QNMs of BSs. Moreover, we extend the results presented in Ref. Macedo et al. (2013b) by computing the QNMs of mini BSs as a function of the central field and compactness. We show that the modes become unstable beyond a critical density corresponding to the maximum mass of the star. In Sec. IV we discuss our results and deliver our final remarks.
Ii Methods
We assume that the perturbation functions can be described by the following set of equations:
(1) 
where is a tortoise radial coordinate, is the frequency of the field, is an dimensional vector representing the perturbations, and is an matrix which can possibly depend on . If the components of the perturbation are decoupled, is a diagonal matrix. In some scenarios, it may be possible to obtain a set of decoupled equations for the perturbations by performing a canonical transformation within a Hamiltonian framework Brizuela (2015). However, the only restriction to here is that it has a diagonal form for , where denotes the characteristic range of the interaction between the components of .
ii.1 Direct integration method
This method consists of integrating the differential equations within two different regions: near the origin and far from the star, with the proper QNM boundary conditions. The problem of finding the QNM frequencies reduces to finding the proper values of for which the solutions obtained integrating from the origin and from infinity are linearly dependent. This method was first used by Chandrasekhar and Detweiler Chandrasekhar and Detweiler (1975), and it was then applied to uniform density stars to compute the least damped (lowest imaginary part) QNMs Chandrasekhar and Ferrari (1991).
The boundary conditions at the origin can be written as^{1}^{1}1The potential usually diverges at the origin due to a centrifugal term, and that is the reason for the term in Eq. (2), where is the angular number of the waves. See, e.g., Ref. Cardoso et al. (2014) and the references therein. We note that, if the central object is a black hole, QNMs are defined by ingoing waves at the event horizon. In this case, a different set of boundary conditions should be imposed in place of Eq. (2).
(2) 
where the are constant dimensional vectors and the upper summation limit is chosen such that the boundary conditions converge to the required accuracy. Substituting Eq. (2) into Eq. (1) and expanding around leads to a recursion relation for the coefficients , such that all of them can be written as functions of the coefficients , which is a collection of independent numbers. In this way, since the system is linear, we can form a set of independent solutions by integrating from the origin by choosing the vector to be, e.g., , and . The general solution can be achieved by a linear combination of the independent solutions, namely,
(3) 
where denotes the th independent solution of Eq. (1), obtained by integrating it from the origin.
The boundary conditions at infinity can be written as
(4) 
where
(5) 
Depending on and on the value of the frequency , we select the value of the sign in Eq. (4) to suit the particular problem. For example, when , the system allows perturbations which are wavelike at infinity, [with being complex]; when , the perturbations are bounded and we require an exponential damping (see, e.g., Table II of Ref. Macedo et al. (2013b)). The are constant dimensional vectors and the upper summation limit in Eq. (4) has to be chosen according to the precision required. Similar to the integration from the origin discussed above, one can construct a set of independent solutions by integrating the perturbation equations from infinity and by choosing different values for . The solution integrated from infinity can be written as the following linear combination:
(6) 
where denotes the th independent solution of Eq. (1). The QNM solutions are such that and are linearly dependent. We see that the QNM frequencies can be found through
(7) 
where . The conditions (7) generate a system of equations for the coefficients and for . Since the system is linear, we can set one of the ’s to unity; say, for instance, . We then use the remaining equations to find the rest of the coefficients as functions of . The remaining equation is then used to find the QNM frequencies.
ii.2 Continuedfraction method
The continuedfraction method is a very powerful technique, with many applications in physics. In the context of QNMs, it was first used by Leaver Leaver (1985), and has been extensively studied by many authors since then Berti et al. (2009). One of the main difficulties in computing QNMs in the frequency domain is the divergence of the wave functions at large distances. The continuedfraction method works extremely well in some cases because it maps the divergent boundary behavior in a specific recurrence relation.
For stellar structures, the continuedfraction method should actually be used in combination with a direct integration. Outside the star (in the vacuum region), one can recast the solution in terms of a continued fraction. However, in order to guarantee that all the proper boundary conditions are satisfied, we have to match continuously this outer solution with the one obtained by integrating the differential equations from the origin directly.
Reference Pani (2013) presents a method for computing the modes through continued fractions, assuming that one can find a series representation in the entire domain. In particular, the method of Ref. Pani (2013) requires that the background spacetime is known analytically. However, in many astrophysical scenarios the metric is constructed only numerically, which makes it difficult to solve the equations through a Frobenius method in the entire domain. Fortunately, even in the presence of interacting fluids and fields, the spacetime might still be described analytically, at least outside an effective radius. This is the case of (slowly rotating) stellar configurations, black holes surrounded by matter, or selfgravitating solitons like BSs. In these cases, it is possible to obtain a continuedfraction representation of the solutions outside an effective radius, and then to match this solution with the one obtained by integrating the equations from requiring regularity at the origin. In this way, we can assure that both boundary conditions are satisfied.
Let us expose the procedure to obtain a series representation outside the effective radius. We again assume that the system of equations can be described by Eq. (1) and that there exists an expansion of the components of the wave vector as follows:
(8) 
where the index denotes the th component of the vector functions and , in which is chosen such that the series solution (8) is convergent Benhar et al. (1999). The vector function is chosen such that the vector satisfies the proper boundary conditions at infinity Konoplya and Zhidenko (2011)—typically, its components are proportional to . Substituting the expansion (8) into the differential equation (1) leads to recurrence relations for the coefficients Leaver (1985); Pani (2013). These could, in principle, be term recurrence relations, which can be reduced to a threeterm recurrence relation by recursive Gaussian elimination steps. The recurrence can be solved, writing all the coefficients —with —in terms of . Bellow, we explain the procedure.
Let us illustrate this by assuming that we end up with a fourterm recurrence relation. We have
(9) 
where , , and are invertible matrices. is the vector whose components are the coefficients in the expansion (8). We can reduce the fourterm recurrence relation to a threeterm one by using a matrixvalued Gaussian elimination step Pani (2013). Using
(10)  
(11)  
(12)  
(13)  
(14) 
one can show, through Eq. (9), that the tilde matrices satisfy the following threeterm recurrence relation:
(15) 
Defining a ladder matrix with the following property,
(16) 
and using Eq. (15), we obtain the following equation:
(17) 
Equation (17) may be solved recursively. We can start at some large value of —say —impose , and then, moving backwards in , determine all ’s. The result has the form of a continued fraction, justifying the name of the method.^{2}^{2}2For instance, if we have just one second order equation, the solution of the recurrence relation (17) takes the form (18) where we used the notation .
Using Eqs. (16) and (17), we obtain all ’s as functions of parameters, given by . The expansion (8), with the determined coefficients, gives the solution to be used outside the star. Therefore, there are independent solutions of Eq. (1), and the general solution is a linear combination of them, similar to Eq. (6). Note that, since the recurrence relations are valid outside the effective radius, the solutions in the interior cannot be described by this procedure.
We still have to impose that the wave function satisfies the proper boundary condition at the origin. For this purpose, we can use the same procedure as in Sec. II.2, obtaining , given by Eq. (3).
The QNM frequencies are found by requiring that the wave functions obtained from inside and outside be linearly dependent. Similarly to the direct integration method, we impose Eq. (7). Note that the difference between the two methods is that we construct the outer solutions through the continuedfraction approximations. For relativistic ordinary stars, the continuedfraction method presented above reduced to the one presented in Refs. Benhar et al. (1999); Nollert (1999).
Iii Application: QNMs of a Boson star
Configurations of Ref. Macedo et al. (2013b)  DI  CF  

MiniBS I, N=1  
MiniBS I, N=2  
MiniBS I, N=3  
MassiveBS I, N=1  
MassiveBS I, N=2  
MassiveBS I, N=3 
Here we extend the calculations of Ref. Macedo et al. (2013b), computing the QNMs of a BS through the improved continuedfraction method. BSs are starlike configurations formed by selfinteracting gravitating complex scalar fields Kaup (1968); Schunck and Mielke (2003). In the case of BSs, outside an effective radius, the perturbations of the spacetime are well described by a decoupled system—two equations describing the polar and axial parts of gravitational perturbations and two equations describing the perturbations of the complex scalar field.^{3}^{3}3As discussed in Ref. Macedo et al. (2013b), although the two equations governing the perturbations for the scalar field are coupled to each other even outside the star, one can always decouple them by a suitable field redefinition Gleiser (1988). The external spacetime is very well approximated by the Schwarzschild geometry and, therefore, the continuedfraction expressions for the wave functions outside the star are the same as in the Schwarzschild case. For the gravitational part of the perturbations, one can construct the continuedfraction relation using the ReggeWheeler equation and use it to obtain the polar and axial gravitational perturbations Leaver (1985); Chandrasekhar (1985). For the scalar field perturbations, we can use the results of Ref. Konoplya and Zhidenko (2005).
For the axial sector of the perturbations, the application of the above method is direct since the equation for the axial perturbations is already in the form of Eq. (1). Moreover, the scalar field perturbations only couple to the polar sector of the perturbations Kojima et al. (1991); Yoshida et al. (1994).
For the polar sector of the perturbations, the method can be applied using the first order differential equations for the gravitational perturbations presented in Ref. Macedo et al. (2013b) in the following way.

We solve the differential equations from the origin, choosing the three free parameters such that we obtain the three independent solutions for the metric and scalar field perturbations.

The vector can be constructed using either the continuedfraction or the direct integration method. Since the gravitational and scalar parts of the perturbations are independent, they can be computed separately;

Finally, with the vector , we find the QNMs by matching it at with , through Eq. (7).
Note that, for the polar case, since we are using the Zerilli function to describe the gravitational perturbations, the point should be such that the background scalar field at that point is negligible.
We can now apply the above procedure to compute the QNMs of a BS. A comparison between results obtained with the direct integration procedure, as in Ref. Macedo et al. (2013b), and with the continuedfraction method for some BS configurations is shown in Table 1. Because of to the small imaginary part of the modes, the divergence of the QNMs is weak and the agreement between the two methods is remarkable. In fact, this good agreement between the two methods is expected because the direct integration method is suitable for computing the least damped modes.
In the case of BSs, the advantage of the continuedfraction method over the direct integration is in the construction of the background solutions. For the direct integration method to work properly, we need to integrate from a point in which has approximately a constant value (or zero). This integration from numerical infinity can be strongly contaminated by numerical errors, such that one generally has to compute the sum in Eq. (4) up to, typically, , and the background solutions need to be constructed up to a point far away from the star. On the other hand, the only requirement for the continuedfraction method to work properly is that the background scalar field is small enough at , and therefore that the continuedfraction solutions can be constructed in a point relatively close to the star.
By using the continuedfraction method, we can now compute the modes of miniBSs as a function of the central density and radius of the star, defined as the radial point which contains of the star’s total mass Macedo et al. (2013b). The results are shown in Figs. 1 and 2. The monopole mode (, Fig. 1) exists due to the coupling between the scalar and gravitational perturbations Yoshida et al. (1994); Macedo et al. (2013b). In Fig. 2, we observe that the quadrupole mode is stable in the considered range of the central density. On the other hand, the monopole mode becomes unstable (; cf. Fig. 3) for central densities , where is the central density corresponding to the maximum mass configuration Schunck and Mielke (2003). This is in accordance with previous works which considered the radial stability of BSs Lee and Pang (1989); Kusmartsev et al. (1991); Hawley and Choptuik (2000).
Iv Discussion and final remarks
We extended two different methods for computing the QNMs of compact objects in the presence of interacting fields. One of the methods is based on a direct integration of the differential equations and the other is based on an analytical continuedfraction scheme for describing the perturbations. The methods can be applied in many astrophysical situations, e.g., to compute the QNMs of relativistic stars surrounded by extra (electromagnetic or exotic) fields.
We applied the direct integration and continuedfraction methods to compute the polar QNMs of BSs. For the modes analyzed, the continuedfraction method is in excellent agreement with the direct integration one. Moreover, we computed the modes as a function of the central value of the scalar field, showing that the monopole mode becomes unstable for configurations beyond the maximum mass.
The techniques presented in this paper can be useful for the computation of QNMs in complicated configurations. One limitation of the continuedfraction method, however, is manifest when dealing with lesser compact interacting fields, i.e., configurations with the longer length scale . In that case, the divergence of the fields usually grows rapidly as increases, and numerical errors can contaminate the accuracy of the solution. One possible way to circumvent that problem is to mix the direct integration method with the continuedfraction one—for instance, one can use the continuedfraction expansion as an outer boundary condition for the integration from infinity, thus decreasing the errors introduced when using the expansion (4).
Acknowledgements.
C.M. and L.C. would like to thank Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), and Fundação Amazônia de Amparo a Estudos e Pesquisas do Pará (FAPESPA). P.P. was supported by the European Community through the IntraEuropean Marie Curie Contract No. AstroGRAphy2013623439 and by FCTPortugal through the Project No. IF/00293/2013. V.C. acknowledges financial support provided under the European Union’s FP7 ERC Starting Grant “The dynamics of black holes: testing the limits of Einstein’s theory” Grant Agreement No. DyBHo–256667, and H2020 ERC Consolidator Grant “Matter and strongfield gravity: New frontiers in Einstein’s theory” Grant Agreement No. MaGRaTh–64659. C.M. thanks the Perimeter Institute for the kind hospitality. Research at the Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development and Innovation. This work was supported by the NRHEP 295189 FP7PEOPLE2011IRSES Grant, and by FCTPortugal through Projects No. PTDC/FIS/098025/2008, No. PTDC/FIS/098032/2008, No. PTDC/FIS/116625/2010, and No. CERN/FP/116341/2010. This work was also partly funded through H2020 ERC Consolidator Grant “Matter and strong field gravity: New frontiers in Einsteinâs theory” Grant Agreement No. MaGRaThâ64659.References
 Berti et al. (2009) E. Berti, V. Cardoso, and A. O. Starinets, Quasinormal modes of black holes and black branes, Class.Quant.Grav. 26, 163001 (2009), arXiv:0905.2975 [grqc] .
 Kokkotas and Schmidt (1999) K. D. Kokkotas and B. G. Schmidt, Quasinormal modes of stars and black holes, Living Rev.Rel. 2, 2 (1999), arXiv:grqc/9909058 [grqc] .
 Kojima (1987) Y. Kojima, Stellar resonant oscillations coupled to gravitational waves, Prog.Theor.Phys. 77, 297–309 (1987).
 Pons et al. (2002) J. A. Pons, E. Berti, L. Gualtieri, G. Miniutti, and V. Ferrari, Gravitational signals emitted by a point mass orbiting a neutron star: Effects of stellar structure, Phys.Rev. D65, 104021 (2002), arXiv:grqc/0111104 [grqc] .
 Macedo et al. (2013a) C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Into the lair: gravitationalwave signatures of dark matter, Astrophys. J. 774, 48 (2013a), arXiv:1302.2646 .
 Barausse et al. (2014) E. Barausse, V. Cardoso, and P. Pani, Can environmental effects spoil precision gravitationalwave astrophysics? Phys. Rev. D 89, 104059 (2014), arXiv:1404.7149 [grqc] .
 Herdeiro and Radu (2014) C. A. R. Herdeiro and E. Radu, Kerr black holes with scalar hair, Phys. Rev. Lett. 112, 221101 (2014), arXiv:1403.2757 [grqc] .
 Herdeiro and Radu (2015) C. A. R. Herdeiro and E. Radu, Asymptotically flat black holes with scalar hair: a review, (2015), arXiv:1504.08209 [grqc] .
 Brito et al. (2015a) R. Brito, V. Cardoso, and P. Pani, Superradiance, Lect. Notes Phys., Vol. 906 (Springer International Publishing, Switzerland, 2015) arXiv:1501.06570 [grqc] .
 Hod (2012) S. Hod, Stationary Scalar Clouds Around Rotating Black Holes, Phys. Rev. D 86, 104026 (2012), arXiv:1211.3202 [grqc] .
 Benone et al. (2014) C. L. Benone, L. C. B. Crispino, C. Herdeiro, and E. Radu, KerrNewman scalar clouds, Phys. Rev. D 90, 104024 (2014), arXiv:1409.1593 [grqc] .
 Barranco et al. (2012) J. Barranco, A. Bernal, J. C. Degollado, A. DiezTejedor, M. Megevand, et al., Schwarzschild black holes can wear scalar wigs, Phys.Rev.Lett. 109, 081102 (2012), arXiv:1207.2153 [grqc] .
 Zilhão et al. (2015) M. Zilhão, H. Witek, and V. Cardoso, Nonlinear interactions between black holes and Proca fields, Class. Quant. Grav. 32, 234003 (2015), arXiv:1505.00797 [grqc] .
 Duncan and Thompson (1992) R. C. Duncan and C. Thompson, Formation of very strongly magnetized neutron stars  implications for gammaray bursts, Astrophys. J. 392, L9 (1992).
 Konno et al. (1999) K. Konno, T. Obata, and Y. Kojima, Deformation of relativistic magnetized stars, Astron. Astrophys. 352, 211–216 (1999), arXiv:grqc/9910038 [grqc] .
 Yazadjiev (2012) S. S. Yazadjiev, Relativistic models of magnetars: Nonperturbative analytical approach, Phys. Rev. D 85, 044030 (2012), arXiv:1111.3536 [grqc] .
 Bertone and Fairbairn (2008) G. Bertone and M. Fairbairn, Compact Stars as Dark Matter Probes, Phys.Rev. D77, 043515 (2008), arXiv:0709.1485 [astroph] .
 Brito et al. (2015b) R. Brito, V. Cardoso, and H. Okawa, Accretion of dark matter by stars, Phys. Rev. Lett. 115, 111301 (2015b), arXiv:1508.04773 [grqc] .
 Brito et al. (2016a) R. Brito, V. Cardoso, C. F. B. Macedo, H. Okawa, and C. Palenzuela, Interaction between bosonic dark matter and stars, Phys. Rev. D 93, 044045 (2016a), arXiv:1512.00466 [astroph.SR] .
 Kaup (1968) D. J. Kaup, KleinGordon Geon, Phys. Rev. 172, 1331–1342 (1968).
 Schunck and Mielke (2003) F. Schunck and E. Mielke, General relativistic boson stars, Class.Quant.Grav. 20, R301–R356 (2003), arXiv:0801.0307 [astroph] .
 Liebling and Palenzuela (2012) S. L. Liebling and C. Palenzuela, Dynamical Boson Stars, Living Rev.Rel. 15, 6 (2012), arXiv:1202.5809 [grqc] .
 Brito et al. (2016b) R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, Proca stars: Gravitating BoseEinstein condensates of massive spin 1 particles, Phys. Lett. B 752, 291–295 (2016b), arXiv:1508.05395 [grqc] .
 Degollado and Herdeiro (2014) J. C. Degollado and C. A. R. Herdeiro, Wiggly tails: a gravitational wave signature of massive fields around black holes, Phys. Rev. D 90, 065019 (2014), arXiv:1408.2589 [grqc] .
 Pani (2013) P. Pani, Advanced Methods in BlackHole Perturbation Theory, Int. J. Mod. Phys. A 28, 1340018 (2013), arXiv:1305.6759 [grqc] .
 Benhar et al. (1999) O. Benhar, E. Berti, and V. Ferrari, The Imprint of the equation of state on the axial w modes of oscillating neutron stars, Mon.Not.Roy.Astron.Soc. 310, 797–803 (1999), arXiv:grqc/9901037 [grqc] .
 Nollert (1999) H.P. Nollert, TOPICAL REVIEW: Quasinormal modes: the characteristic ‘sound’ of black holes and neutron stars, Class.Quant.Grav. 16, R159–R216 (1999).
 Macedo et al. (2013b) C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Astrophysical signatures of boson stars: quasinormal modes and inspiral resonances, Phys. Rev. D 88, 064046 (2013b), arXiv:1307.4812 [grqc] .
 Brizuela (2015) D. Brizuela, A generalization of the Zerilli master variable for a dynamical spherical spacetime, Class. Quant. Grav. 32, 135005 (2015), arXiv:1505.05278 [grqc] .
 Chandrasekhar and Detweiler (1975) S. Chandrasekhar and S. L. Detweiler, The quasinormal modes of the Schwarzschild black hole, Proc.Roy.Soc.Lond. A344, 441–452 (1975).
 Chandrasekhar and Ferrari (1991) S. Chandrasekhar and V. Ferrari, On the NonRadial Oscillations of a Star III. A Reconsideration of the Axial Modes, Proceedings: Mathematical and Physical Sciences 434, pp. 449–457 (1991).
 Cardoso et al. (2014) V. Cardoso, L. C. B. Crispino, C. F. B. Macedo, H. Okawa, and P. Pani, Light rings as observational evidence for event horizons: longlived modes, ergoregions and nonlinear instabilities of ultracompact objects, Phys. Rev. D 90, 044069 (2014), arXiv:1406.5510 [grqc] .
 Leaver (1985) E. Leaver, An Analytic representation for the quasi normal modes of Kerr black holes, Proc. Roy. Soc. Lond. A 402, 285–298 (1985).
 Konoplya and Zhidenko (2011) R. A. Konoplya and A. Zhidenko, Quasinormal modes of black holes: From astrophysics to string theory, Rev. Mod. Phys. 83, 793–836 (2011).
 Gleiser (1988) M. Gleiser, Stability of Boson Stars, Phys.Rev. D38, 2376 (1988).
 Chandrasekhar (1985) S. Chandrasekhar, The mathematical theory of black holes (Oxford University Press, Oxford, 1985).
 Konoplya and Zhidenko (2005) R. Konoplya and A. Zhidenko, Decay of massive scalar field in a Schwarzschild background, Phys. Lett. B 609, 377–384 (2005), arXiv:grqc/0411059 [grqc] .
 Kojima et al. (1991) Y. Kojima, S. Yoshida, and T. Futamase, Nonradial pulsation of a boson star. 1: Formulation, Prog.Theor.Phys. 86, 401–410 (1991).
 Yoshida et al. (1994) S. Yoshida, Y. Eriguchi, and T. Futamase, Quasinormal modes of boson stars, Phys.Rev. D50, 6235–6246 (1994).
 Lee and Pang (1989) T. Lee and Y. Pang, Stability of Mini  Boson Stars, Nucl.Phys. B315, 477 (1989).
 Kusmartsev et al. (1991) F. V. Kusmartsev, E. W. Mielke, and F. E. Schunck, Gravitational stability of boson stars, Phys.Rev. D43, 3895–3901 (1991), arXiv:0810.0696 [astroph] .
 Hawley and Choptuik (2000) S. H. Hawley and M. W. Choptuik, Boson stars driven to the brink of black hole formation, Phys. Rev. D 62, 104024 (2000), arXiv:grqc/0007039 [grqc] .