QNMs of slowly rotating Einstein–Bumblebee black hole

We have studied the quasinormal modes (QNMs) of a slowly rotating black hole with Lorentz-violating parameter in Einstein–Bumblebee gravity. We analyse the slow rotation approximation of the rotating black hole in the Einstein–Bumblebee gravity, and obtain the master equations for scalar perturbation, vector perturbation and axial gravitational perturbation, respectively. Using the matrix method and the continuous fraction method, we numerically calculate the QNM frequencies. In particular, for scalar field, it shows that the QNMs up to the second order of rotation parameter have higher accuracy. The numerical results show that, for both scalar and vector fields, the Lorentz-violating parameter has a significant effect on the imaginary part of the QNM frequencies, while having a relatively smaller impact on the real part of the QNM frequencies. But for axial gravitational perturbation, the effect of increasing the Lorentz-violating parameter ℓ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell $$\end{document} is similar to that of increasing the rotation parameter a~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{a}}$$\end{document}.


I. INTRODUCTION
Lorentz symmetry is a very basic and important spacetime symmetry, which is not only the basis of quantum field theory and the standard model of modern particle physics, but also plays a fundamental role in relativity.Motivated by establishing the quantum gravity, and the evidence from high energy cosmic rays, the possibility of Lorentz symmetry breaking (LSB) on Planck energy scale is expected and has been extensively discussed [1].The suggestions for sources of Lorentz violation (LV) included string field theory [2,3], loop quantum gravity theory [4], noncommutative field theories [5,6], quantum gravity inspired spacetime foam scenarios [7], and varying couplings [8,9].However, in the low energy scale, this Lorentz symmetry breaking and its effects should be largely supressed, otherwise it would be inconsistent with many current ground-based experimental observations.Thus, over the past few decades, the search for quantum gravity effects in the low energy scale attached much attention.Especially, the standard-model extension is a candidate which offers a broad theoretical foundation for testing Lorentz symmetry.
The Lorentz symmetry is broken when a vector field has a nonzero vacuum expectation value.A straightforward effective theory of gravity with a standard model extension LV is the Bumblebee model, where the Lorentz violation results from the dynamics of a single vector or axial-vector field B µ , also known as the bumblebee field.Kostelecky and Samuel first used the bumblebee gravitational model to investigate the effects of spontaneous LV [10,11].Bertolami and Páramos derived the vacuum solutions of this gravity model, including the purely radial LSB, the radial/temporal LSB or the axial/temporal LSB [12,13].In recent years, Casana et al. presented a Schwarzschild-like bumblebee black hole solution [14].After that, a series of other spherically symmetric solutions within the framework of the bumblebee gravity theory have been discovered, including a traversable wormhole solution [17], the solution with the global monopole [18], the cosmological constant [19], or the Einstein-Gauss-Bonnet term [20].A Kerr-like exact solution in Bumblebee gravity was obtained by Ding et al. [21], which does not seem to satisfy the corresponding field equation [22].However, we find that this solution can be satisfied under some certain conditions.Linear perturbation of black holes play a important role in physics.Many astrophysical processes can be regarded as small derivations from the black hole background spacetime.For example, QNMs is considered as a good description of the late stage of the merge of binary black holes or the gravitational collapse.In order to calculate the QNMs, one need obtain the decoupled second-order differential equation in the frequency domain.For spherically symmetric spacetimes, one can always decompose the perturbed metric into spherical harmonics.And this process generally ensures that the angular part and the radial part can be automatically decoupled, while the odd-parity perturbation and the evenparity perturbation are also automatically separated [23][24][25][26].But this decomposition is no longer valid in rotating spacetime.One important work is Teukolsky discovered that a class of perturbation equations of Kerr spacetime can be decoupled in Newman-Penrose formalism [27], the underlying reason is that the Kerr spacetime is of type-D in the Petrov classification.However, it is very difficult to decouple the perturbation equations in other rotating spacetimes, such as Kerr-Newman metric or metric for most modified gravity theories.For Proca perturbation of Kerr spacetime, the separability of perturbation equation also seems not easy [28].In recently, Frolov et.al.successfully solved this problem [29].
Recently, a new idea was provided by Pani et al., that is, if the slowly rotating backgrounds are close enough to a spherical symmetric spacetime, then the separability of the perturbation equations becomes possible [30][31][32][33][34][35].This method originated from Kojima's work on the perturbations of slowly rotating neutron stars [36].Under the slowly rotating limit, considering that the dimensionless rotation parameter ã = a/M of the background spacetime is sufficiently small, using the harmonic basis to expand the metric perturbation, the field equations can always reduce to a coupled ordinary differential equations.Different from the spherically symmetric case, in the limit of slowly rotating there exist a "selection rule", i.e. the perturbations with a given parity and multipolar index l are coupled to three parts [30]: the perturbations with opposite parity and index l±1 at first order of ã, the perturbations with the same parity and the same index l up to second order of ã2 , and the perturbations with the same parity and index l ± 2 at the second order of ã.
The aim of this manuscript is to apply the slow rotation method to Einstein-bumblebee gravity theories, and obtain the massive scalar perturbation, the Proca field perturbation and the axial gravitational perturbation for rotating Einstein-bumblebee black hole, respectively.We also calculate the QNMs of these perturbations by numerical methods, and study how the LV parameter affect the QNMs.The manuscript is organized as follows.In Sec.II, we briefly review the Einstein-Bumblebee theory and the slow rotation technique.In Sec.III, we derive the Schrödinger-like equation for the massive scalar perturbation, the proca perturbation and the axial gravitational perturbation for slowly rotating Einstein-Bumblebee black hole, respectively.In Sec.IV, using both the matrix method and the continued fraction method, we calculated the QNM frequencies of rotating Einstein-Bumblebee black holes.Sec.V is devoted to summary and discussion.

A. The Field Equations
Here we briefly review the Einstein-Bumblebee gravity model and the black hole solutions under this gravity.This model is known as an example that extends the standard formalism of general relativity.Under a suitable potential, the bumblebee vector field B µ acquires a nonzero vacuum expectation value and induces a spontaneous Lorentz symmetry breaking in the gravitational sector.The action for a single bumblebee field B µ coupled to gravity can be described as [14][15][16] ̺ is the real coupling constant which controls the nonminimal gravity interaction to bumblebee field B a , and the bumblebee field strength is defined by The potential V is chosen to provide a nonvanishing vacuum expectation value for bumblebee field B a , and has a minimum at B a B a ± b 2 , where b is a real positive constant.Hence the potential can be expressed as [14] and this potential can drive a nonzero vaccum value Taking the variation of Eq. ( 1) yields the vacuum gravitational equation and the equation of motion for the bumblebee field where κ = 8πG N , T B ab is the bumblebee energy momentum tensor, and V ′ denotes ∂V (y)/∂y at y = B a B a ± b 2 .Using the trace of Eq. ( 4), we obtain Now further assume that the bumblebee field is fixed to be B a = b a , then the particular form of the potential is irrelevant and and the gravitational field Eq. ( 6) is equivalent to In the following subsection, the satisfied of Eq. (10) determines that whether the metric is the exact solution of the vacuum Einstein-Bumblebee action.
B. Slowly Rotating Einstein-Bumblebee Black Hole Recently, using the condition b a b a = constant, and then considering the bumblebee field b a = (0, b 0 ρ √ ∆ , 0, 0), Ding et al. find the Kerr-like solution for Einstein-Bumblebee theory can be written as [21] where and ℓ = ̺b 2 0 is the Lorentz-violating parameter.This metric is the solution for rotating spacetime with a radial bumblebee field.It is obvious that the metric becomes the usual Kerr metric when ℓ → 0 and the static Einstein-Bumblebee metric when ã → 0. Unfortunately, the metric (11) seems not the exact-solution [22].Considering that ã is a sufficiently small parameter, and expanding the metric to the second order of rotation parameter ã, we obtain Note that ∆ is equivalent to ∆ in the second order of ã, and which can determine the singular points r + and r − as the roots of ∆ ∆ = (r − r + )(r − r − ). ( where We find that this metric yields some components of the field equation do not satisfied.However, we find the non vanishing terms for field equations are all proportional to ℓ 2 ã2 .Under certain conditions, this metric can be considered as an approximate solution of the Einstein-Bumblebee Eq. (10).The detailed discussions are presented in Appendix A.

III. PERTURBATIONS OF SLOWLY ROTATION EINSTEIN-BUMBLEBEE BLACK HOLE
In this section, using the slow rotation approximation, we analysis the perturbations of Einstein-Bumblebee black hole.For massive scalar field perturbation, our analysis up to the second order of ã.And for massive vector field perturbation or odd-parity gravitational perturbation, our analysis only up to the first order of ã.
Here we briefly introduce the notations of the spherical harmonic basis [37].First we define two unnormalized and orthogonal co-vectors v and n as v a = (−1, 0, 0, 0), n a = (0, 1, 0, 0), (18) the projection operator onto the sphere surface and the spatial Levi-Civita tensor, ǫ abc ≡ v d ǫ dabc , where ǫ trθφ = r 2 sin θ. using the scalar spherical harmonics Y lm = Y lm (θ, ϕ), the pure-spin vector harmonics are given by and the pure-spin tensor harmonics are given by These spherical harmonic basis will be used to decompose vector or tensor perturbations.For axially symmetric background, the perturbation with different values of m are decoupled, therefore we ignore the index m in the following discussions.And the linearized perturbed equations always imply a sum over (l, m).Some identities will be used, which are given by [36] cos where A. Massive scalar perturbation Considering that the scalar field coupling to the bumblebee field is neglected, then the massive Klein-Gordon equation reads where m s = µ represents the scalar field's mass.We decompose the field in spherical harmonics: Since in axially symmetric background, perturbations with different values of m are decoupled, the index of m can be neglected.Substituting the metric (15) into Eq.( 27) and expanding the equation up to the second order of ã, we obtain We find that the parameters A l and D l are constructed by Ψ l and its derivatives, and present their specific expression in Appendix B. It shows that in Eq. ( 29) the radial and the angular sections are still coupled together, hence one needs to decouple the angular section to obtain a purely radial equation.Expanding Eq. ( 29) up to the first order of ã, one can obtain where dr/dx = F = 1−2M/r, and the effective potential is to represent the potential expanded up to the first order of ã.When the Lorentz-violating parameter ℓ vanishes, the effective potential will be consistent with Eq. ( 23) in Ref. [30].Now we separate the angular part of Eq. ( 29) up to the second order of ã.Using the identity Eq. ( 24) and taking the inner product of the Eq. ( 29) with the conjugate of scalar spherical harmonics on the sphere, we can get the parameters D l±2 is made up of Ψ l±2 and its derivatives, which are written as Note that the explicit expression of W l±2 is unnecessary.The reason is that all D l±2 terms is proportional to ã2 , one should only consider D l±2 at the zeroth order of ã.This implies that the functions Ψ l±2 in Eq. ( 33) can be only considered as the solutions of where And the expression of D l±2 at the zeroth order can be obtained, Then, using the explicit form of the coefficients A l±2 and D l±2 , the field equation( 32) has the form as where dr/dr * = f , and the potential is given by This potential coincides with Pani's result when ℓ → 0.
To obtain the single variable equation for Ψ l , we need define where c l reads as Then, the wave-like equation for scalar field perturbation is given by and which can be analyzed by standard methods.

B. Massive Vector Perturbation
The field equations of a massive vector field, or known as Proca field, is where A a is the vector potential, µ = m v / is the mass of the vector field and . We also consider this vector field satisfied the Lorentz condition ∇ a A a = 0.With the help of the vector spherical harmonic basis (20), we expand the perturbed vector δA a as, where Λ = l(l + 1), and P l , R l , S l , Q l are all scalar functions of coordinate (t, r).The functions P l , R l , S l belong to the polar sector and the function Q l belongs to the axial sector.Moreover, assume an harmonic timedependence in time, for instance, The perturbation of the four components of Eq.( 42) and the Lorentz condition can be naturally separated into three groups: The index I = {t, r, L} represents the t component, the r component of Eq. ( 45), and the expansion of the Lorentz condition ∇ a A a = 0, respectively.The explicit expression of the coefficients in Eqs.(45) to (47) are presented in Appendix C.
For Eq. ( 45), the angular part can be separated by taking the inner product to Y * lm .But for Eqs. ( 46) and (47), in order to completely separate the angular part, we need construct a new vector Πa as Πa ≡ (0, 0, δΠ θ , δΠ ϕ ) . (48) Taking the inner product with Πa and also using the identities Eqs. ( 22)-( 25), then the radial equations can be obtained.We list these radial equations in Appendix C, see Eqs. (C24)-(C26).
First we consider the static Proca field equations.
where the operator D2 is introduced as It shows that in this case, two sectors are naturally decoupled.Using Eq. ( 52) to eliminate R l in Eq. ( 51), one can obtain which is a fourth order differential equation.Note that we can classify the eigenvectors of the system according to the three degrees of freedom of the vector A a , i.e., the three polarizations, which include one scalar type polarization and two vector type polarizations.And the electric mode of the vector potential A a has one scalar type polarization and one vector type polarization.The reason of the above equation as a fourth order differential equation is that, S l contains both scalar type polarization and vector type polarization.Now we consider the situation of expanding the field equations to the first order of ã.The polar sector of the field equations gives where S po 1 and S po 2 are the source terms, and the explicit form of which are presented in Appendix C. While the axial sector of the field equations gives where S ax is the source term and given in Appendix C. It shows that the expressions of S po 1 , S po 2 and S ax are proportional to ã, therefore it is only necessary to consider the components of S, such as Q l or R l , to the zeroth order of ã.
Note that for the first order rotation approximation, the massive scalar perturbation and the massive vector perturbation can be uniformly written as where s = 0 for scalar perturbations and s = ±1 for vector perturbations with axial parity.However, from the nextsubsection, it shows that in Einstein-Bumblebee theory, the gravitational perturbation with axial parity can not be uniformly written in the form of the above equation.

C. Gravitational Perturbation
We consider that the metric g ab can be decomposed into the background g (0) ab and the perturbation h ab Using the scalar spherical harmonic, the pure-spin vector harmonics (20) and the tensor harmonics (21), one can obtain ten orthogonal spherical harmonic basis [37].And the metric perturbation can be decomposed as [38] where the coefficients A l to K l are scalar functions.
Adopting the well known Regge-Wheeler gauge [23,[38][39][40][41][42][43][44], we set After separating the angular components of Eq. ( 10), ten pure radial equations can be obtained, which are naturally separated into axial parity and polar parity.The two independent components for the axial parity equation are given by Rrϕ : Together with Eqs. ( 64) and ( 63), define U l as then the modified master equation describing the axial gravitational perturbation up to the first order of ã can be written as where From the expression of the above potential, we find that the Lorentz-violating parameter coupled with rotation parameter ã, which implies that the axial parity perturbation equation for the static black hole solution in Einstein-Bumblebee theory is the same as Schwarzschild black hole.

IV. THE EIGENVALUE PROBLEM FOR QUASINORMAL MODES A. Boundary conditions
In this manuscript, we are concerned about how the Lorentz-violating parameter ℓ affects the QNMs.Hence we only investigate the QNMs at the massless limit.At the first order of ã, the relation (17) shows that there is only one horizon at r + = 2M .For scalar field or vector field, the generic wave function Ψ l has the asymptotic behavior as where x is the tortoise coordinate given in Eq. ( 30), l and σ + are defined as But for the gravitational field Eq. ( 67), the generic wave function U l has the asymptotic behavior as Using these asymptotic solutions, in order to apply the numerical method, we can impose Ψ l and U l satisfied the relation as If we consider the case that up to the second order of rotation parameter ã, two horizons of the rotation black hole are determined by the Eq. ( 17).At this case the asymptotic behavior of Z l can be written as where According to this asymptotic behavior, Z l can be assumed to be written as In the following sections, we will analyze QNMs by using the matrix method and the continued fraction method, respectively.

B. Matrix method for quasinormal modes
To compute the QNMs, here we briefly describe the matrix method presented by Lin et.al.[45][46][47][48].By taking into account the boundary conditions mentioned in the previous subsection, we perform a coordinate transformation so the region of QNMs calculation becomes y ∈ [0, 1].Considering the boundary condition, assuming the wave function can be reconstructed as then the boundary condition at the event horizon and the spatial infinity becomes This boundary condition ensures that the resulting matrix equation is homogenous, as will be seen below.It can be proved that all perturbation field equations can be rewritten as where the functions C(y, ω) can be derived by substituting the behavior of the wave function into the corresponding field equations.For instance, together with Eqs. ( 76) and (77), one can substitute Eq. ( 75) into ( 41) and obtain the equation in the form of Eq. ( 79).Note that all C j (j = 0, 1, 2) are linear functions of ω, so C j can be decomposed as C j (y, ω) = C j,0 (y) + ωC j,1 (y).
Using the matrix method to discretize Eq. ( 79), we introduce equally spaced grid points into the internal [0, 1].By expanding the function χ(y) around each grid point using the Taylor series, the corresponding differential matrices can be constructed.Thus, Eq. ( 79) is rewritten as an algebraic equation in matrix form where M 0 and M 1 are matrices consisting of the functions C j and the corresponding differential matrices.Calculating these matrices [45,46], and then the solve of QNMs becomes a simple algebraic solution problem.

C. Continued fraction method
Since the groundbreaking work by Leaver [49], the continued fraction method is an accurate method in determining the QNMs.In this method, the eigenfunction can be expressed as a series whose coefficients adhere to a finite term recurrence relation.In order to obtain a more concise recurrence relation, in this subsection we set M = 1/2.
A solution to the perturbation equation that expands at the event horizon can be written in the following form: For the Eq.(59) controlling the scalar and electromagnetic fields, the expansion coefficients are defined by a three-term recursion relations The recurrence coefficients α n , β n , and γ n are simple functions consisting of n and other differential equation parameters, the explicit forms are as follows − 4 l[ lΛ + lmã(i + 3 lω)] + 16 lω(i + 2 lω), (85) The quasi-normal frequencies can be obtained by solving the algebraic Eq. ( 83) at sufficiently large n for any initial value d 0 .Without loss of generality, we set d 0 = 1.
For the case of gravitational perturbations Eq. ( 67), the recursion relations appear to be more complicated.The explicit form of these relations are given by where Ci,j , as well as α n , β n , γ n , σ n , τ n and δ n , are functions consisting of n and other parameters.The explicit form are presented in Appendix D.

D. Numerical results
Using the matrix method and the continued fraction method, we numerically calculated QNMs and show the result in this subsection.In matrix method, we set N = 18, to ensure that the relative error becomes smaller than 10 −5 .In continuous fraction method, we computed the 60th order for scalar or electromagnetic perturbations and 16th order for gravitational perturbation.The reason is that the latter equation is more complex and consumes too much computing resources.However, the numerical results shows that the difference between two methods is smaller than 10 −4 .In the calculating of QNMs, we set M = 1.In this paper we only provide data for the l = m = 2 modes, which expected to be the most astrophysically relevant.The detailed calculation data can be found in Appendix E.

Static of Bumblebee modes
First, we consider the static black hole in the Einstein-Bumblebee theory.When ã → 0, the metric (11) becomes the static Schwarzschild-like black hole solution given by Ref. [14].Note that this static Einstein-Bumblebee solution is an exactly solution without any approximation.From the determinant of the metric (11), in order to maintain the Lorentz signature, the parameter ℓ should satisfied ℓ > −1 [16].For the gravitational perturbation of the static Einstein-Bumblebee black hole, the perturbation equation is the same as that in Schwarzschild case, hence the QNM frequencies of the gravitational field are completely independent of the Lorentz-violating parameter ℓ.We calculated the QNMs for massless scalar perturbation, the electromagnetic perturbation with the range 0 < ℓ < 1, respectively.Fig. 1-2 show the QNM frequencies for the scalar field and the electromagnetic field in the complex plane as calculated via matrix method.

A consistency check for massless scalar perturbation
Before the calculation of the QNMs of the rotation Einstein-Bumblebee black hole, we consider the special case of scalar field at ℓ = 0, i.e., the Kerr black hole.The QNMs of the Kerr black holes are widely discussed [45,[49][50][51].We compute the QNMs by slow rotation approximation at first and second order, and then compare the results with the exact Kerr results provided by [51].Fig. 3 show that, the QNM frequencies results calculated via slow rotation approximation matches the exact massless scalar perturbation of a Kerr black hole, for l = m = 2 mode.The figures reveal that, including the real and imaginary parts, the difference between our approach and the exact values is about 1% up to ã = a/M ∼ 0.4.It is obvious that the first order approximation for the QNM frequencies begins to fail at smaller values of ã, and the second order approximation has much higher accuracy.And it also shows that the frequencies of QNMs do not deviate significantly from the standard result until ã ∼ 0.2.

QNMs of Slowly rotating Einstein-bumblebee black hole
Now we turn to the slowly rotating Einstein-Bumblebee case with non-vanishing ã and ℓ.The rotation parameter ã is limited in the range 0 ≤ ã ≤ 0.3.
For massless scalar field, we computed the QNM frequencies under both the first order and the second order approximations of the rotation.Fig. 4 show that, for determined ℓ, the imaginary part of the second-order approximation shows a more obvious change with the increase of ã.While for determined ã, the absolution of the imaginary part will decrease with the increase of ℓ.
For the electromagnetic and gravitational fields, we calculate the QNM frequencies under the first order approximation of the rotation parameter ã.Figs. 5 and  6 show that, in the rotating Einstein-Bumblebee black hole, with the increase of ℓ of the Lorentz-violating parameter, the real part of the QNM frequencies will in- crease, but the absolution of the imaginary part will decrease.
Fig. 6 also shows quite intuitively that with a relatively small rotation parameter, ã < 0.1, the Lorentzviolating parameter can hardly have a significant effect on the QNM frequencies.The reason is that for gravitational fields, the Lorentz-violating parameter only affects the first-order rotational correction term of the effective potential.Analyzing the effective potential of the gravitational field shows that following the definition the form of the gravitational field master Eq. (67) will remain the same as in the Kerr case [34].It should be noted that this equivalence is a coincidence at the first order of ã.Whether this coincidence would be satisfied at the higher order needs further exploration.
The results of QNMs for gravitational perturbations also show that, the parameters ã and ℓ affecting the QNMs are degenerate.This implies that it is hard to determine whether the deviation of QNMs is coming from the rotation or the LV, and the effect of the bumblebee field may be misinterpreted as being due to the rotation.

V. CONCLUSIONS AND EXTENSIONS
In this paper, we investigate the calculations of QNMs of the rotation Einstein-Bumblebee black hole, which is an approximation solution of stationary and axisymmetric black hole of Einstein-Bumblebee theory.In this solution, the bumblebee field only has a pure radial vacuum energy expectation and assumed to be as b µ = (0, b r , 0, 0).The strength of the LV can be determined by a coupling constant ℓ.For the rotation solution, expanding the metric to the second order of rotation parameter ã, we obtain the first and second order differential equations for the massive scalar perturbation.And we also give the Proca field perturbation equation and the odd parity of gravitational perturbation equation at the first order of ã, respectively.It shows that the first order scalar field equation and the first order Proca field equation can be written in a form that includes these two perturbations.
Using the matrix method and the continued fraction method, we calculated the QNMs for different fields.Our results show that, the QNMs errors obtained by the two methods are negligible.It clearly show that, as we expected, compare with the exactly result in Kerr solution, the second order approximation has a higher accuracy than the first order approximation.For scalar, vector and gravitational perturbation, we find that the real part of the QNM frequencies are not sensitive to the change of Lorentz-violating parameter ℓ.However, the results  shows that with the increase of ℓ, the absolute values of the imaginary part of the QNM frequencies decrease, which indicate that the perturbation will decay slower.
It should be noted that the rotation metric we consider corresponding to the bumblebee field has only the radial component.In this theory the Lorentz-violating parameter ℓ is determined by both coupling constant ̺ and non-zero VEV b a .Actually, the bumblebee field can be assumed to be b µ = (0, b(r), b(θ), 0) [52] or b µ = (b(t), b(r), 0, 0) [53], and the corresponding metrics can be obtained.In these cases, how LV affects the QNMs is also a open problem.
Another point worth noting is that the isospectrality of axial and polar perturbations has been shown in many spacetimes, for example, the Schwarzschild metric and the Reissner Nordström metric [54].The numerical results of Pani et al. provide a strong evidence that the electromagnetic and the gravitational modes of slowly rotation Kerr metric are isospectral at the first order of ã [32].However, the isospectrality is violated in some modified theories of gravity, the violation of isospectrality, such as Lovelock gravity [55], Chern-Simons gravity [56] and loop quantum gravity [57].For Einstein-bumblebee gravity, since we do not know how to construct the master equation for even parity of electromagnetic or gravitational perturbations, whether the bumblebee black hole has the isospectral property is still an open question.
Appendix C: Coefficients and source terms of the Proca perturbation In this appendix we present the explicit expressions of the parameters and the relations for the massive vector perturbation.The explicit form of the coefficients appearing in Eqs. ( 45)-( 47) is B Appendix E: QNM FREQUENCY TABLES

FIG. 1 :
FIG. 1: Complex massless scalar frequencies for the n = 0 mode for varying values of ℓ.

FIG. 2 :
FIG. 2: Complex massless electromagnetic frequencies for the n = 0 mode for varying values of ℓ.

FIG. 3 :
FIG. 3: Comparison between the exact Kerr result and the results obtained by slow rotation at first or second order for varying values of a with l = m = 2 and M = 1.

20 FIG. 4 :
FIG.4: Real and imaginary parts of the l = m = 2 massless scalar mode for varying values of a, where the dashed line represents the first-order approximation of the rotation parameter and the solid line represents the second-order approximation of the rotation parameter.

2 FIG. 5 :
FIG. 5: Real and imaginary parts of the l = m = 2 massless electromagnetic mode for varying values of a/M that belong to the axial sector.

FIG. 6 :
FIG. 6: Real and imaginary parts of the l = m = 2 gravitational mode for varying values of a/M .

TABLE I :
Comparison of the n = 0, l = m = 2 mode massless scalar QNM frequencies calculated by the matrix method and the continued fraction method under the first order slow rotation approximation.

TABLE II :
Comparison of the n = 0, l = m = 2 mode massless electromagnetic QNM frequencies calculated by the matrix method and the continued fraction method under the first order slow rotation approximation.

TABLE III :
Comparison of the n = 0, l = m = 2 mode gravitational QNM frequencies calculated by the matrix method and the continued fraction method under the first order slow rotation approximation.