QNMs of slowly rotating Einstein–Bumblebee black hole

Abstract

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 $\ell$ is similar to that of increasing the rotation parameter ${\tilde{a}}$.

1 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_\mu$, 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 even-parity 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 ${\tilde{a}}=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\pm 1$ at first order of ${\tilde{a}}$, the perturbations with the same parity and the same index l up to second order of ${\tilde{a}}^2$, and the perturbations with the same parity and index $l\pm 2$ at the second order of ${\tilde{a}}$.

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 Sect. 2, we briefly review the Einstein–Bumblebee theory and the slow rotation technique. In Sect. 3, 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 Sect. 4, using both the matrix method and the continued fraction method, we calculated the QNM frequencies of rotating Einstein–Bumblebee black holes. Section 5 is devoted to summary and discussion.

2 The theoretical framework

2.1 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_\mu$ 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_\mu$ coupled to gravity can be described as [14,15,16]

$$\begin{aligned} {\mathcal {S}}&=\int d^4x \sqrt{-g}\left[ \frac{1}{16\pi G_N}\left( {\mathcal {R}}+\varrho B^aB^b {\mathcal {R}}_{ab} \right) \right. \nonumber \\&\quad \left. -\frac{1}{4}B^{ab}B_{ab}-V \right] +{\mathcal {L}}_M, \end{aligned}$$

(1)

$\varrho$ is the real coupling constant which controls the non-minimal gravity interaction to bumblebee field $B_a$, and the bumblebee field strength is defined by

$$\begin{aligned} B_{ab}=\partial _a B_b-\partial _b B_a. \end{aligned}$$

(2)

The potential V is chosen to provide a nonvanishing vacuum expectation value for bumblebee field $B_a$, and has a minimum at $B_aB^a\pm b^2$, where b is a real positive constant. Hence the potential can be expressed as [14]

$$\begin{aligned} V=V(B_aB^a\pm b^2), \end{aligned}$$

(3)

and this potential can drive a nonzero vaccum value $\langle B^a\rangle =b^a$ with $b_ab^a=\mp b^2$.

Taking the variation of Eq. (1) yields the vacuum gravitational equation and the equation of motion for the bumblebee field

$$\begin{aligned}&{\mathcal {R}}_{ab}-\frac{1}{2}g_{ab}{\mathcal {R}}=\kappa T^{B}_{ab}, \end{aligned}$$

(4)

$$\begin{aligned}&\nabla ^aB_{ab}=2V'B_b-\frac{\varrho }{\kappa }B^aR_{ab}. \end{aligned}$$

(5)

where $\kappa =8\pi G_N$, $T^{B}_{ab}$ is the bumblebee energy momentum tensor,

$$\begin{aligned} T^{B}_{ab}&= B_{ac}B^c_b-\frac{1}{4}g_{ab}B^{cd}B_{cd}-g_{ab}V+2B_aB_bV' \nonumber \\&\quad \times \frac{\varrho }{\kappa }\left[ \frac{1}{2}g_{ab}B^{c}B^{d}R_{cd}-B_aB^cR_{cb}-B_bB^cR_{ca}\right. \nonumber \\&\quad \left. +\frac{1}{2}\nabla _c\nabla _a\left( B^cB_b\right) +\frac{1}{2}\nabla _c\nabla _b\left( B^cB_a\right) \right. \nonumber \\&\quad \left. -\frac{1}{2}\nabla ^2\left( B_aB_b\right) -\frac{1}{2}g_{ab}\nabla _c\nabla _d\left( B^cB^d\right) \right] , \end{aligned}$$

(6)

and $V'$ denotes $\partial V(y)/\partial y$ at $y=B^aB_a\pm b^2$. Using the trace of Eq. (4), we obtain

$$\begin{aligned} {\mathcal {R}}_{ab}&=\kappa T^{B}_{ab}+2\kappa g_{ab}V-\kappa g_{ab}B^cB_cV'\nonumber \\&\quad +\frac{\varrho }{4}g_{ab}\nabla ^2\left( B^cB_c \right) +\frac{\varrho }{2}g_{ab}\nabla _c\nabla _d\left( B^cB^d \right) . \end{aligned}$$

(7)

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 $V=V'=0$. Define

$$\begin{aligned} {\bar{R}}_{ab}&={\mathcal {R}}_{ab}-\kappa b_{ac}b^c_b+\frac{\kappa }{4}g_{ab}b^{cd}b_{cd}+\varrho b_ab^c{\mathcal {R}}_{cb}\nonumber \\&\quad +\varrho b_bb^c{\mathcal {R}}_{ca}-\frac{\varrho }{2}g_{ab}b^cb^d{\mathcal {R}}_{cd}+\bar{{\mathcal {B}}}_{ab}, \end{aligned}$$

(8)

$$\begin{aligned} \bar{{\mathcal {B}}}_{ab}&=\frac{\varrho }{2}\left[ \nabla ^2\left( b_ab_b\right) -\nabla _c\nabla _a\left( b^cb_b\right) -\nabla _c\nabla _b\left( b^cb_a\right) \right] , \end{aligned}$$

(9)

and the gravitational field Eq. (6) is equivalent to

$$\begin{aligned} {\bar{R}}_{ab}=0. \end{aligned}$$

(10)

In the following subsection, the satisfied of Eq. (10) determines that whether the metric is the exact solution of the vacuum Einstein–Bumblebee action.

2.2 Slowly rotating Einstein–Bumblebee black hole

Recently, using the condition $b^ab_a=constant$, and then considering the bumblebee field $b_a=(0,b_0\frac{\rho }{\sqrt{\Delta }},0,0)$, Ding et al. find the Kerr-like solution for Einstein–Bumblebee theory can be written as [21]

$$\begin{aligned} ds^2&=-\left( 1-\frac{2M r}{\rho ^2}\right) dt^2-\frac{4M r a\sqrt{1+\ell }\sin ^2\theta }{\rho ^2}dtd\varphi \nonumber \\&\quad +\frac{\rho ^2}{\Delta }dr^2+\rho ^2d\theta ^2+\frac{{\mathcal {A}}\sin ^2\theta }{\rho ^2}d\varphi ^2 \end{aligned}$$

(11)

where

$$\begin{aligned} \rho ^2&=r^2+(1+\ell )a^2\cos ^2\theta , \end{aligned}$$

(12)

$$\begin{aligned} \Delta&=\frac{r^2-2Mr}{1+\ell }+a^2, \end{aligned}$$

(13)

$$\begin{aligned} {\mathcal {A}}&=\left[ r^2+(1+\ell )a^2 \right] ^2-\Delta (1+\ell )^2a^2\sin ^2\theta . \end{aligned}$$

(14)

and $\ell =\varrho b_0^2$ 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 $\ell \rightarrow 0$ and the static Einstein–Bumblebee metric when ${\tilde{a}}\rightarrow 0$. Unfortunately, the metric (11) seems not the exact-solution [22]. Considering that ${\tilde{a}}$ is a sufficiently small parameter, and expanding the metric to the second order of rotation parameter ${\tilde{a}}$, we obtain

$$\begin{aligned} ds^2&=g_{ab}^{(2)}dx^adx^b \nonumber \\&= -\left( 1-\frac{2M}{r}+\frac{2a^2M(1+\ell )\cos ^2\theta }{ r^3}\right) dt^2 \nonumber \\&\quad -\frac{4aM\sqrt{1+\ell }\sin ^2\theta }{r}dtd\varphi +\frac{(1+\ell )\rho ^2}{\tilde{\Delta }}dr^2\nonumber \\&\quad +\left[ r^2+a^2(1+\ell )\cos ^2\theta \right] d\theta ^2 \nonumber \\&\quad +\sin ^2\theta \left[ r^2+a^2\left( 1+\ell \right) \left( 1+\frac{2M}{r}\sin ^2\theta \right) \right] d \varphi ^2, \end{aligned}$$

(15)

Note that $\tilde{\Delta }$ is equivalent to $\Delta$ in the second order of ${\tilde{a}}$, and which can determine the singular points $r_{+}$ and $r_{-}$ as the roots of $\tilde{\Delta }$

$$\begin{aligned} \tilde{\Delta }=(r-r_{+})(r-r_{-}). \end{aligned}$$

(16)

where

$$\begin{aligned} r_{+}=2M-(1+\ell )\frac{a^2}{2M}, \quad r_{-}=(1+\ell )\frac{a^2}{2M}. \end{aligned}$$

(17)

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 $\ell ^2{\tilde{a}}^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.

3 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 ${\tilde{a}}$. And for massive vector field perturbation or odd-parity gravitational perturbation, our analysis only up to the first order of ${\tilde{a}}$. 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

$$\begin{aligned} v_a=(-1,0,0,0), \quad n_a=(0,1,0,0), \end{aligned}$$

(18)

the projection operator onto the sphere surface

$$\begin{aligned} \Omega _{ab}=r^2 \text {diag}(0,0,1,\sin ^2\theta ), \end{aligned}$$

(19)

and the spatial Levi–Civita tensor, $\epsilon _{abc}\equiv v^d\epsilon _{dabc}$, where $\epsilon _{tr\theta \phi }=r^2\sin \theta$. using the scalar spherical harmonics $Y^{lm}=Y^{lm}(\theta ,\varphi )$, the pure-spin vector harmonics are given by

$$\begin{aligned} Y_a^{E,lm}&= r\nabla _a Y^{lm}, \nonumber \\ Y_a^{B,lm}&= r\epsilon _{ab}^c n^b\nabla _c Y^{lm}, \nonumber \\ Y_a^{R,lm}&= n_aY^{lm}, \end{aligned}$$

(20)

and the pure-spin tensor harmonics are given by

$$\begin{aligned} T_{ab}^{T0,lm}&= \Omega _{ab}Y^{lm}, \nonumber \\ T_{ab}^{L0,lm}&= n_an_b Y^{lm}, \nonumber \\ T_{ab}^{E1,lm}&= rn_{(a}\nabla _{b)}Y^{lm}, \nonumber \\ T_{ab}^{B1,lm}&= rn_{(a}\epsilon _{b)c}^d n^c\nabla _d Y^{lm}, \nonumber \\ T_{ab}^{E2,lm}&= r^2\left( \Omega _a^c\Omega _b^d-\frac{1}{2}\Omega _{ab}\Omega ^{cd}\right) \nabla _c\nabla _dY^{lm}, \nonumber \\ T_{ab}^{B2,lm}&= r^2\Omega _{(a}^c\epsilon _{b)e}^d n^e\nabla _c\nabla _d Y^{lm}. \end{aligned}$$

(21)

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]

$$\begin{aligned} \cos \theta Y^{l}&={\mathcal {Q}}_{l+1}Y^{l+1}+{\mathcal {Q}}_{l}Y^{l-1}, \end{aligned}$$

(22)

$$\begin{aligned} \sin \theta \partial _\theta Y^{l}&={\mathcal {Q}}_{l+1}lY^{l+1}-{\mathcal {Q}}_{l}(l+1)Y^{l-1}, \end{aligned}$$

(23)

$$\begin{aligned} \cos ^2\theta Y^{l}&=\left( {\mathcal {Q}}^2_{l+1}+{\mathcal {Q}}^2_{l}\right) Y^{l} \nonumber \\&\quad +{\mathcal {Q}}_{l+1}{\mathcal {Q}}_{l+2}Y^{l+2} +{\mathcal {Q}}_{l}{\mathcal {Q}}_{l-1}Y^{l-2}, \end{aligned}$$

(24)

$$\begin{aligned} \cos \theta \sin \theta \partial _\theta Y^{l}&=\left[ l{\mathcal {Q}}^2_{l+1}-(l+1){\mathcal {Q}}^2_{l}\right] Y^{l}\nonumber \\&\quad +{\mathcal {Q}}_{l+1}{\mathcal {Q}}_{l+2}lY^{l+2}\nonumber \\&\quad -{\mathcal {Q}}_{l}{\mathcal {Q}}_{l-1}(l+1)Y^{l-2} , \end{aligned}$$

(25)

where

$$\begin{aligned} {\mathcal {Q}}_{l}=\sqrt{\frac{l^2-m^2}{4l^2-1}}. \end{aligned}$$

(26)

3.1 Massive scalar perturbation

Considering that the scalar field coupling to the bumblebee field is neglected, then the massive Klein-Gordon equation reads

$$\begin{aligned}&\frac{1}{\sqrt{-g}}\partial _a\left( \sqrt{-g}g^{ab}\partial _b\phi \right) =\mu ^2\phi , \end{aligned}$$

(27)

where $m_s=\mu \hbar$ represents the scalar field’s mass. We decompose the field in spherical harmonics:

$$\begin{aligned} \phi =\sum _{lm}\frac{\Psi _{l}(r)}{\sqrt{r^2+a^2}}e^{-i\omega t}Y^{lm}(\theta ,\varphi ). \end{aligned}$$

(28)

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 ${\tilde{a}}$, we obtain

$$\begin{aligned} A_lY^{l}+a^2 D_l\cos ^2\theta \cdot Y^{l}=0, \end{aligned}$$

(29)

We find that the parameters $A_l$ and $D_l$ are constructed by $\Psi _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 ${\tilde{a}}$, one can obtain

$$\begin{aligned} \frac{d^2}{dx^2}\Psi _l^{(1)}+V^{(1)}_l\Psi _l^{(1)}=0, \end{aligned}$$

(30)

where $dr/dx={\mathcal {F}}=1-2M/r$, and the effective potential is

$$\begin{aligned} V^{(1)}_l&=(1+\ell )\left( \omega ^2-\sqrt{1+\ell }\frac{4aMm\omega }{r^3}\right) \nonumber \\&\quad -{\mathcal {F}}\left[ \frac{2M}{r^3}+(1+\ell )\left( \frac{l(l+1)}{r^2}+\mu ^2\right) \right] , \end{aligned}$$

(31)

to represent the potential expanded up to the first order of ${\tilde{a}}$. When the Lorentz-violating parameter $\ell$ 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 ${\tilde{a}}$. 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

$$\begin{aligned}&A_l+a^2\left[ \left( {\mathcal {Q}}^2_{l+1}+{\mathcal {Q}}^2_l\right) D_l+{\mathcal {Q}}_{l-1}{\mathcal {Q}}_lD_{l-2}\right. \nonumber \\&\quad \left. +{\mathcal {Q}}_{l+2}{\mathcal {Q}}_{l+1}D_{l+2}\right] =0, \end{aligned}$$

(32)

the parameters $D_{l\pm 2}$ is made up of $\Psi _{l\pm 2}$ and its derivatives, which are written as

$$\begin{aligned} D_{l\pm 2}=\frac{1}{r^2(2M-r)}\left( \frac{d^2}{d x^2}\Psi _{l\pm 2}+W_{l\pm 2}\Psi _{l\pm 2} \right) . \end{aligned}$$

(33)

Note that the explicit expression of $W_{l\pm 2}$ is unnecessary. The reason is that all $D_{l\pm 2}$ terms is proportional to ${\tilde{a}}^2$, one should only consider $D_{l\pm 2}$ at the zeroth order of ${\tilde{a}}$. This implies that the functions $\Psi _{l\pm 2}$ in Eq. (33) can be only considered as the solutions of

$$\begin{aligned} \frac{d^2}{dx^2}\Psi _{l\pm 2}^{(0)}+V^{(0)}_{l\pm 2}\Psi _{l\pm 2}^{(0)}=0. \end{aligned}$$

(34)

where $V^{(0)}=V^{(1)}|_{a=0}$. And the expression of $D_{l\pm 2}$ at the zeroth order can be obtained,

$$\begin{aligned} D_{l\pm 2}^{(0)}=\frac{(1+\ell )(\omega ^2-\mu ^2)}{r^3}\Psi _{l\pm 2}^{(0)}, \end{aligned}$$

(35)

Then, using the explicit form of the coefficients $A_{l\pm 2}$ and $D_{l\pm 2}$, the field equation (32) has the form as

$$\begin{aligned}&\frac{d^2}{dr_*^2}\Psi _l^{(2)}+V_l^{(2)}\Psi _l^{(2)}=\frac{a^2{\mathcal {F}}(1+\ell )^2(\mu ^2-\omega ^2)}{r^2}\nonumber \\&\quad \times \left[ {\mathcal {Q}}_{l+1}{\mathcal {Q}}_{l+2}\Psi _{l+2}^{(0)}+{\mathcal {Q}}_{l-1}{\mathcal {Q}}_{l}\Psi _{l-2}^{(0)} \right] , \end{aligned}$$

(36)

where $dr/dr_*=f$,

$$\begin{aligned} f=\frac{(1+\mathcal {\ell })\tilde{\Delta }}{r^2+a^2}, \end{aligned}$$

(37)

and the potential is given by

$$\begin{aligned} V_l^{(2)}&=V^{(1)}_l+\frac{a^2}{r^6}\Big [-24M^2+2Mr\big (6-2l(l+1)(1+\ell ) \nonumber \\&\quad -3\ell -2r^2(1+\ell ) \mu ^2+r^2(1+\ell )^2\omega ^2 \big ) \nonumber \\&\quad +r^2\left( l-1+2\ell -l\ell ^2+m^2(1+\ell )^2-l^2(\ell ^2-1)\right. \nonumber \\&\quad \left. +r^2\mu ^2-r^2\left( \omega ^2+\ell ^2\mu ^2-\ell ^2\omega ^2 \right) \right) \nonumber \\&\quad -r^4{\mathcal {F}}(1+\ell )^2(\mu ^2-\omega ^2)\left( {\mathcal {Q}}_{l}^2+{\mathcal {Q}}_{l+1}^2\right) \Big ]. \end{aligned}$$

(38)

This potential coincides with Pani’s result when $\ell \rightarrow 0$. To obtain the single variable equation for $\Psi _l$, we need define

$$\begin{aligned} Z_l=\Psi _l^{(2)}+a^2c_l \Psi _{l-2}^{(2)}-a^2c_{l+2}\Psi _{l+2}^{(2)}, \end{aligned}$$

(39)

where $c_l$ reads as

$$\begin{aligned} c_l=\frac{(1+\ell )}{2(2l-1)}\left( \mu ^2-\omega ^2\right) {\mathcal {Q}}_{l-1}{\mathcal {Q}}_{l}. \end{aligned}$$

(40)

Then, the wave-like equation for scalar field perturbation is given by

$$\begin{aligned} \frac{d^2}{dr_*^2}Z_l+V_l^{(2)}Z_l=0, \end{aligned}$$

(41)

and which can be analyzed by standard methods.

3.2 Massive vector perturbation

The field equations of a massive vector field, or known as Proca field, is

$$\begin{aligned} \Pi ^b=\nabla _aF^{ab}-\mu ^2A^b=0, \end{aligned}$$

(42)

where $A^a$ is the vector potential, $\mu =m_v/\hbar$ is the mass of the vector field and $F_{ab}=2\nabla _{[a}A_{b]}$. We also consider this vector field satisfied the Lorentz condition $\nabla _aA^a=0$. With the help of the vector spherical harmonic basis (20), we expand the perturbed vector $\delta A_a$ as,

$$\begin{aligned} \delta A_a&=-\frac{\textrm{P}_l}{r}v_aY^{lm}+\frac{\textrm{R}_l}{rf}Y^{R,lm}_a\nonumber \\&\quad +\frac{\textrm{S}_l}{r\Lambda }Y^{E,lm}_a -\frac{\sqrt{1+\ell } \textrm{Q}_l}{r\Lambda }Y^{B,lm}_a, \end{aligned}$$

(43)

where $\Lambda =l(l+1)$, and $\textrm{P}_l$, $\textrm{R}_l$, $\textrm{S}_l$, $\textrm{Q}_l$ are all scalar functions of coordinate (t, r). The functions $\textrm{P}_l$, $\textrm{R}_l$, $\textrm{S}_l$ belong to the polar sector and the function $\textrm{Q}_l$ belongs to the axial sector. Moreover, assume an harmonic time-dependence in time, for instance,

$$\begin{aligned} \textrm{P}_l(t,r)=e^{-i \omega t}\textrm{P}_l(r). \end{aligned}$$

(44)

The perturbation of the four components of Eq. (42) and the Lorentz condition can be naturally separated into three groups:

$$\begin{aligned} \delta \Pi _{I}&=\left( A^{(I)}_l+{\tilde{A}}^{(I)}_l\cos \theta +D^{(I)}_l\cos ^2\theta \right) Y^l\nonumber \\&\quad +\left( B^{(I)}_l+{\tilde{B}}^{(I)}_l \cos ^2\theta \right) \sin \theta \frac{\partial }{\partial \theta }Y^l=0, \end{aligned}$$

(45)

$$\begin{aligned} \delta \Pi _\theta&=\left( \alpha _l+\rho _l\sin ^2\theta \right) \frac{\partial }{\partial \theta }Y^l-i m \beta _l \frac{Y^l}{\sin \theta }\nonumber \\&\quad +\left( \eta _l+\sigma _l\cos \theta \right) \sin \theta Y^l=0, \end{aligned}$$

(46)

$$\begin{aligned} \frac{\delta \Pi _\varphi }{\sin \theta }&=\left( \beta _l+\gamma _l\sin ^2\theta \right) \frac{\partial }{\partial \theta }Y^l+i m \alpha _l\frac{Y^l}{\sin \theta }\nonumber \\&\quad +\left( \zeta _l+\lambda _l\cos \theta \right) \sin \theta Y^l=0. \end{aligned}$$

(47)

The index $I=\{t,r,L\}$ represents the t component, the r component of Eq. (45), and the expansion of the Lorentz condition $\nabla _aA^a=0$, respectively. The explicit expression of the coefficients in Eqs. (45)–(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 $\tilde{\Pi }_a$ as

$$\begin{aligned} \tilde{\Pi }_a\equiv \left( 0,0,\delta \Pi _\theta ,\delta \Pi _\varphi \right) . \end{aligned}$$

(48)

Taking the inner product with $Y^{a*}_{E,l'm'}$ and $Y^{a*}_{B,l'm'}$, respectively,

$$\begin{aligned}&\int \tilde{\Pi }_a Y^{a*}_{E,l'm'}d\Omega , \end{aligned}$$

(49)

$$\begin{aligned}&\int \tilde{\Pi }_a Y^{a*}_{B,l'm'}d\Omega , \end{aligned}$$

(50)

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. When ${\tilde{a}}\rightarrow 0$, Eqs. (C24)–(C26) yield

$$\begin{aligned}&\hat{{\mathcal {D}}}_2 \textrm{R}_l-\frac{2{\mathcal {F}}}{r^2}\left( 1-\frac{3M}{r}\right) \left[ \textrm{R}_l-(1+\ell )\textrm{S}_l \right] =0, \end{aligned}$$

(51)

$$\begin{aligned}&\hat{{\mathcal {D}}}_2 \textrm{S}_l+\frac{2 \Lambda {\mathcal {F}}}{r^2}\textrm{R}_l=0, \end{aligned}$$

(52)

$$\begin{aligned}&\hat{{\mathcal {D}}}_2 \textrm{Q}_l=0, \end{aligned}$$

(53)

where the operator $\hat{{\mathcal {D}}}_2$ is introduced as

$$\begin{aligned} \hat{{\mathcal {D}}}_2=\frac{d^2}{dx^2}+(1+\ell )\left[ \omega ^2-{\mathcal {F}}\left( \frac{\Lambda }{r^2}+\mu ^2\right) \right] . \end{aligned}$$

(54)

It shows that in this case, two sectors are naturally decoupled. Using Eq. (52) to eliminate $\textrm{R}_l$ in Eq. (51), one can obtain

$$\begin{aligned} \hat{{\mathcal {D}}}_2\left[ \frac{r^2}{{\mathcal {F}}}\hat{{\mathcal {D}}}_2\textrm{S}_l \right] \!-\!\left( 1-\frac{3M}{r}\right) \left[ 2\hat{{\mathcal {D}}}_2 +(1\!+\!\ell )\frac{4\Lambda {\mathcal {F}}}{r^2} \right] \textrm{S}_l=0, \end{aligned}$$

(55)

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, $\textrm{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 ${\tilde{a}}$. The polar sector of the field equations gives

$$\begin{aligned}&\hat{{\mathcal {D}}}_2 \textrm{R}_l-\frac{2{\mathcal {F}}}{r^2}\left( 1-\frac{3M}{r}\right) \bigg [\textrm{R}_l-(1+\ell )\textrm{S}_l \bigg ]={\mathcal {S}}^{\text {po}}_1, \end{aligned}$$

(56)

$$\begin{aligned}&\hat{{\mathcal {D}}}_2\textrm{S}_l+\frac{2\Lambda {\mathcal {F}}}{r^2}\textrm{R}_l={\mathcal {S}}^{\text {po}}_2, \end{aligned}$$

(57)

where ${\mathcal {S}}^\text {po}_1$ and ${\mathcal {S}}^\text {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

$$\begin{aligned} \hat{{\mathcal {D}}}_2 \textrm{Q}_l-(1+\ell )^{\frac{3}{2}}\frac{4a Mm\omega }{r^3}\textrm{Q}_l={\mathcal {S}}^{\text {ax}}. \end{aligned}$$

(58)

where ${\mathcal {S}}^\text {ax}$ is the source term and given in Appendix C. It shows that the expressions of ${\mathcal {S}}^\text {po}_1$, ${\mathcal {S}}^\text {po}_2$ and ${\mathcal {S}}^\text {ax}$ are proportional to ${\tilde{a}}$, therefore it is only necessary to consider the components of ${\mathcal {S}}$, such as $\textrm{Q}_l$ or $\textrm{R}_l$, to the zeroth order of ${\tilde{a}}$.

Note that for the first order rotation approximation, the massive scalar perturbation and the massive vector perturbation can be uniformly written as

$$\begin{aligned} \hat{{\mathcal {D}}}_2 \Psi _l-\frac{2M}{r^3}\left[ (1-s^2){\mathcal {F}}+2a m(1+\ell )^\frac{3}{2}\omega \right] \Psi _l=0. \end{aligned}$$

(59)

where $s=0$ for scalar perturbations and $s=\pm 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.

3.3 Gravitational perturbation

We consider that the metric $g_{ab}$ can be decomposed into the background $g^{(0)}_{ab}$ and the perturbation $h_{ab}$

$$\begin{aligned} g_{ab}=g^{(0)}_{ab}+h_{ab}. \end{aligned}$$

(60)

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]

$$\begin{aligned} h_{ab}^{lm}&= \textrm{A}_l~v_av_bY^{lm}+2\textrm{B}_l~v_{(a}Y_{b)}^{E,lm}+2\textrm{C}_l~v_{(a}Y_{b)}^{B,lm}\nonumber \\&\quad +2\textrm{D}_l~v_{(a}Y_{b)}^{R,lm}+\textrm{E}_l~T^{T0,lm}_{ab} +\textrm{F}_l~T_{ab}^{E2,lm}\nonumber \\&\quad +\textrm{G}_l~T^{B2,lm}_{ab}+2\textrm{H}_l~T^{E1,lm}_{ab}+2\textrm{J}_l~T^{B1,lm}_{ab}\nonumber \\&\quad +\textrm{K}_l~T^{L0,lm}_{ab} . \end{aligned}$$

(61)

where the coefficients $\textrm{A}_l$ to $\textrm{K}_l$ are scalar functions. Adopting the well known Regge-Wheeler gauge [23, 38,39,40,41,42,43,44], we set

$$\begin{aligned} \textrm{B}_l=\textrm{F}_l=\textrm{H}_l=\textrm{G}_l=0, \end{aligned}$$

(62)

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

$$\begin{aligned}&{\bar{R}}_{t\varphi }: \frac{\Lambda {\mathcal {F}}}{2} \left( 2\textrm{C}_l'+r\textrm{C}_l''-3i\omega \textrm{J}_l -ir\omega \textrm{J}_l'\right) +\frac{aMm\sqrt{1+\ell }}{r^2}\nonumber \\&\quad \times \left[ \frac{2-\Lambda }{{\mathcal {F}}}\omega \textrm{C}_l +i\left( \frac{12M}{r^2}\textrm{J}_l-\frac{4+\Lambda }{r}\textrm{J}_l+2{\mathcal {F}}\textrm{J}_l'\right) \right] \nonumber \\&\quad +\frac{\Lambda (4M-r\Lambda )}{2r^2}\textrm{C}_l=0, \end{aligned}$$

(63)

$$\begin{aligned}&{\bar{R}}_{r\varphi }: \frac{(1+\ell )\Lambda }{2{\mathcal {F}}}\left( \frac{2(2-\Lambda )M}{r^2}+\frac{\Lambda -2-r^2\omega ^2}{r} \right) \textrm{J}_l\nonumber \\&\quad +\frac{i\omega (1+\ell )\Lambda }{2{\mathcal {F}}}\left( \textrm{C}_l-r\textrm{C}_l' \right) +\frac{aMm(1+\ell )^{\frac{3}{2}}}{{\mathcal {F}}}\nonumber \\&\quad \times \left[ i\frac{6-\Lambda }{r^3}\textrm{C}_l+\frac{\Lambda }{r^2}\left( 2\omega \textrm{J}_l+i\textrm{C}_l' \right) \right] =0. \end{aligned}$$

(64)

Together with Eqs. (64) and (63), define ${\mathcal {U}}_l$ as

$$\begin{aligned} \textrm{J}_l=\frac{1}{{\mathcal {F}}}\left( 1-2\sqrt{1+\ell }\frac{aMm}{r^3\omega } \right) {\mathcal {U}}_l, \end{aligned}$$

(65)

then the modified master equation describing the axial gravitational perturbation up to the first order of ${\tilde{a}}$ can be written as

$$\begin{aligned} \frac{d^2}{d x^2}{\mathcal {U}}_l+\left[ \omega ^2-{\mathcal {V}}_l\right] {\mathcal {U}}_l=0, \end{aligned}$$

(66)

where

$$\begin{aligned} {\mathcal {V}}_l&={\mathcal {F}}\left( \frac{\Lambda }{r^2}-\frac{6M}{r^3} \right) +\sqrt{1+\ell }\frac{4aMm}{r^3}\nonumber \\&\quad \times \left[ \omega +6{\mathcal {F}}\left( \frac{3r-7M}{r^3\Lambda \omega } \right) \right] . \end{aligned}$$

(67)

From the expression of the above potential, we find that the Lorentz-violating parameter coupled with rotation parameter ${\tilde{a}}$, 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.

4 The eigenvalue problem for quasinormal modes

4.1 Boundary conditions

In this manuscript, we are concerned about how the Lorentz-violating parameter $\ell$ affects the QNMs. Hence we only investigate the QNMs at the massless limit. At the first order of ${\tilde{a}}$, the relation (17) shows that there is only one horizon at $r_{+}=2M$. For scalar field or vector field, the generic wave function $\Psi _l$ has the asymptotic behavior as

$$\begin{aligned} \Psi _l\sim {\left\{ \begin{array}{ll} \left( r-2M\right) ^{-i \tilde{\ell }\sigma _+} &{}\text {for}~~r\rightarrow ~2M ,\\ e^{i\tilde{\ell }\omega x} &{}\text {for}~~r\rightarrow ~+\infty , \end{array}\right. } \end{aligned}$$

(68)

where x is the tortoise coordinate given in Eq. (30), $\tilde{\ell }$ and $\sigma _{+}$ are defined as

$$\begin{aligned} \tilde{\ell }=\sqrt{1+\ell }, \quad \sigma _{+}=2M\omega -\tilde{\ell }\frac{m a}{2M}. \end{aligned}$$

(69)

But for the gravitational field Eq. (67), the generic wave function ${\mathcal {U}}_l$ has the asymptotic behavior as

$$\begin{aligned} {\mathcal {U}}_l\sim {\left\{ \begin{array}{ll} \left( r-2M\right) ^{-i \sigma _+} &{}\text {for}~~r\rightarrow ~2M ,\\ e^{i\omega x} &{}\text {for}~~r\rightarrow ~+\infty . \end{array}\right. } \end{aligned}$$

(70)

Using these asymptotic solutions, in order to apply the numerical method, we can impose $\Psi _l$ and ${\mathcal {U}}_l$ satisfied the relation as

$$\begin{aligned} \Psi _l&= e^{i\tilde{\ell }\omega r}r^{i\tilde{\ell }\left( 2M \omega +\sigma _+\right) }\left( r-2M\right) ^{-i\tilde{\ell }\sigma _+}\psi _l, \end{aligned}$$

(71)

$$\begin{aligned} {\mathcal {U}}_l&= e^{i\omega r}r^{i\left( 2M \omega +\sigma _+\right) }\left( r-2M\right) ^{-i\sigma _+}\psi _l. \end{aligned}$$

(72)

If we consider the case that up to the second order of rotation parameter ${\tilde{a}}$, 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

$$\begin{aligned} Z_l\sim {\left\{ \begin{array}{ll} \left( r-r_+\right) ^{-i \tilde{\ell }\Omega } &{}\text {for}~~r\rightarrow ~r_+ ,\\ e^{i\tilde{\ell }\omega r_*} &{}\text {for}~~r\rightarrow ~+\infty , \end{array}\right. } \end{aligned}$$

(73)

where

$$\begin{aligned} r_{+}=2M-(1+\ell )\frac{a^2}{2M}, \quad \Omega =\left( 4M-r_+\right) \omega -\tilde{\ell } \frac{ma}{2M}. \end{aligned}$$

(74)

According to this asymptotic behavior, $Z_l$ can be assumed to be written as

$$\begin{aligned} Z_l=&e^{i\tilde{\ell }\omega r}r^{i\tilde{\ell }\left( 2M \omega +\Omega \right) }\left( r-r_+\right) ^{-i\tilde{\ell }\Omega }\psi _l. \end{aligned}$$

(75)

In the following sections, we will analyze QNMs by using the matrix method and the continued fraction method, respectively.

4.2 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

$$\begin{aligned} y=1-\frac{r_+}{r}, \end{aligned}$$

(76)

so the region of QNMs calculation becomes $y\in [0,1]$. Considering the boundary condition, assuming the wave function can be reconstructed as

$$\begin{aligned} \chi (y)=y(1-y)\psi _l(y), \end{aligned}$$

(77)

then the boundary condition at the event horizon and the spatial infinity becomes

$$\begin{aligned} \chi (0)=\chi (1)=0. \end{aligned}$$

(78)

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

$$\begin{aligned} {\mathcal {C}}_2(y,\omega )\chi ''(y)+{\mathcal {C}}_1(y,\omega )\chi '(y)+{\mathcal {C}}_0(y,\omega )\chi (y)=0, \end{aligned}$$

(79)

where the functions ${\mathcal {C}}(y,\omega )$ 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 ${\mathcal {C}}_j(j=0,1,2)$ are linear functions of $\omega$, so ${\mathcal {C}}_j$ can be decomposed as ${\mathcal {C}}_j(y,\omega )={\mathcal {C}}_{j,0}(y)+\omega {\mathcal {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 $\chi (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

$$\begin{aligned} \left( {\mathcal {M}}_0+\omega {\mathcal {M}}_1\right) \chi (y)=0, \end{aligned}$$

(80)

where ${\mathcal {M}}_0$ and ${\mathcal {M}}_1$ are matrices consisting of the functions ${\mathcal {C}}_j$ and the corresponding differential matrices. Calculating these matrices [45, 46], and then the solve of QNMs becomes a simple algebraic solution problem.

4.3 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:

$$\begin{aligned} \psi _l=\sum _{n=0}^{\infty }d_n\left( \frac{r-r_+}{r-r_-} \right) ^n. \end{aligned}$$

(81)

For the Eq. (59) controlling the scalar and electromagnetic fields, the expansion coefficients are defined by a three-term recursion relations

$$\begin{aligned}&[1+i(1+\ell )m{\tilde{a}}-2i\tilde{\ell }\omega ]d_1+\big [s^2-1-(1+\ell )(\Lambda -8\omega ^2)\nonumber \\&\quad +4i\tilde{\ell }\omega -(1+\ell )m{\tilde{a}}(i+3\tilde{\ell }\omega )\big ]d_0=0, \end{aligned}$$

(82)

$$\begin{aligned}&d_{n+1}\alpha _n+d_n\beta _n+d_{n-1}\gamma _n=0,n=1,2,3\cdots \end{aligned}$$

(83)

The recurrence coefficients $\alpha _n$, $\beta _n$, and $\gamma _n$ are simple functions consisting of n and other differential equation parameters, the explicit forms are as follows

$$\begin{aligned} \alpha _n&=4(1+n)(1+n+i(1+\ell )m{\tilde{a}}-2i\tilde{\ell }\omega ), \end{aligned}$$

(84)

$$\begin{aligned} \beta _n&=4(s^2-1-2n^2)-8n[1+i(1+\ell )m{\tilde{a}}-4i\tilde{\ell }\omega ]\nonumber \\&\quad -4\tilde{\ell }[\tilde{\ell }\Lambda +\tilde{\ell }m{\tilde{a}}(i+3\tilde{\ell }\omega )]+16\tilde{\ell }\omega (i+2\tilde{\ell }\omega ), \end{aligned}$$

(85)

$$\begin{aligned} \gamma _n&=4(n^2-s^2)+4in[(1+\ell )m{\tilde{a}}-4\tilde{\ell }\omega ]\nonumber \\&\quad +8(1+\ell )(\tilde{\ell }m{\tilde{a}}-2\omega )\omega . \end{aligned}$$

(86)

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

$$\begin{aligned}&d_1=\tilde{{\mathcal {C}}}_{1,0}~d_0,\nonumber \\&d_2 =\tilde{{\mathcal {C}}}_{2,0}~d_0+\tilde{{\mathcal {C}}}_{2,1}~d_1,\nonumber \\&d_3 =\tilde{{\mathcal {C}}}_{3,0}~d_0+\tilde{{\mathcal {C}}}_{3,1}~d_1+\tilde{{\mathcal {C}}}_{3,2}~d_2,\nonumber \\&d_4 =\tilde{{\mathcal {C}}}_{4,0}~d_0+\tilde{{\mathcal {C}}}_{4,1}~d_1+\tilde{{\mathcal {C}}}_{4,2}~d_2+\tilde{{\mathcal {C}}}_{4,3}~d_3, \end{aligned}$$

(87)

$$\begin{aligned}&d_{n+1}\alpha _n+d_{n}\beta _n+d_{n-1}\gamma _n \nonumber \\&\quad +d_{n-2}\sigma _n+d_{n-3}\tau _n+d_{n-4}\delta _n=0, \quad n=4,5,6\cdots \end{aligned}$$

(88)

where $\tilde{{\mathcal {C}}}_{i,j}$, as well as $\alpha _n$, $\beta _n$, $\gamma _n$, $\sigma _n$, $\tau _n$ and $\delta _n$, are functions consisting of n and other parameters. The explicit form are presented in Appendix D.

4.4 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.

4.4.1 Static of Bumblebee modes

First, we consider the static black hole in the Einstein–Bumblebee theory. When ${\tilde{a}}\rightarrow 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 $\ell$ should satisfied $\ell >-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 $\ell$. We calculated the QNMs for massless scalar perturbation, the electromagnetic perturbation with the range $0<\ell <1$, respectively. Figures 1 and 2 show the QNM frequencies for the scalar field and the electromagnetic field in the complex plane as calculated via matrix method.

Fig. 1

Complex massless scalar frequencies for the $n=0$ mode for varying values of $\ell$

Fig. 2

Complex massless electromagnetic frequencies for the $n=0$ mode for varying values of $\ell$

4.4.2 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 $\ell =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

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$

Figure 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 ${\tilde{a}}=a/M\sim 0.4$. It is obvious that the first order approximation for the QNM frequencies begins to fail at smaller values of ${\tilde{a}}$, 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 ${\tilde{a}}\sim 0.2$.

4.4.3 QNMs of slowly rotating Einstein–Bumblebee black hole

Now we turn to the slowly rotating Einstein–Bumblebee case with non-vanishing ${\tilde{a}}$ and $\ell$. The rotation parameter ${\tilde{a}}$ is limited in the range $0\le {\tilde{a}}\le 0.3$.

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

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

Real and imaginary parts of the $l=m=2$ gravitational mode for varying values of a/M

For massless scalar field, we computed the QNM frequencies under both the first order and the second order approximations of the rotation. Figure 4 show that, for determined $\ell$, the imaginary part of the second-order approximation shows a more obvious change with the increase of ${\tilde{a}}$. While for determined ${\tilde{a}}$, the absolution of the imaginary part will decrease with the increase of $\ell$.

For the electromagnetic and gravitational fields, we calculate the QNM frequencies under the first order approximation of the rotation parameter ${\tilde{a}}$. Figures 5 and 6 show that, in the rotating Einstein–Bumblebee black hole, with the increase of $\ell$ of the Lorentz-violating parameter, the real part of the QNM frequencies will increase, but the absolution of the imaginary part will decrease.

Figure 6 also shows quite intuitively that with a relatively small rotation parameter, ${\tilde{a}}<0.1$, the Lorentz-violating 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

$$\begin{aligned} {\bar{a}}=\sqrt{1+\ell }a, \end{aligned}$$

(89)

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 ${\tilde{a}}$. 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 ${\tilde{a}}$ and $\ell$ 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.

5 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_{\mu }=(0,b_r,0,0)$. The strength of the LV can be determined by a coupling constant $\ell$. For the rotation solution, expanding the metric to the second order of rotation parameter ${\tilde{a}}$, 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 ${\tilde{a}}$, 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 $\ell$. However, the results shows that with the increase of $\ell$, 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 $\ell$ is determined by both coupling constant $\varrho$ and non-zero VEV $b^a$. Actually, the bumblebee field can be assumed to be $b_{\mu }=(0,b(r),b(\theta ),0)$ [52] or $b_{\mu }=(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 ${\tilde{a}}$ [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.

Appendices

Appendix A: Conditions for solutions up to second order of ${\tilde{a}}$

In this appendix, we discuss that the metric (15) is the slowly approximation solution of the Einstein–Bumblebee theory, i.e., satisfied the condition Eq. (10).

By resolving the Einstein–Bumblebee gravitational field equations, Ding et. al. presented an Kerr-like solution. But Ref. [22] pointed out that their solution is not the exact solution. However, with sufficiently small Lorentz-violating parameter $\ell$, considering the slow rotation approximation and the parameter ${\tilde{a}}$ is small, we can get an approximation solution. Substituting the metric (15) into the Eq. (10), we find the zero components include

$$\begin{aligned} {\bar{R}}_{tr}={\bar{R}}_{t\theta }={\bar{R}}_{t\varphi }={\bar{R}}_{r\theta }={\bar{R}}_{r\varphi }={\bar{R}}_{\theta \varphi }=0, \end{aligned}$$

(A1)

and the non-zero components are

$$\begin{aligned} {\bar{R}}_{tt}&=\frac{r-2M}{r}\delta , \quad {\bar{R}}_{rr}=\frac{r(1+\ell )}{r-2M}\delta ,\nonumber \\ {\bar{R}}_{\theta \theta }&=-r^2\delta ,\quad {\bar{R}}_{\varphi \varphi }=-r^2\sin ^2\theta ~\delta , \end{aligned}$$

(A2)

where $\delta$ is

$$\begin{aligned} \delta =\frac{a^2\ell ^2}{4r^4}\left( 1+3\cos 2\theta \right) . \end{aligned}$$

(A3)

If $\ell$ and ${\tilde{a}}$ are small parameters, we consider that the non vanishing ${\bar{R}}_{tt}$ and ${\bar{R}}_{\theta \theta }$ can be neglected. Hence the metric (15) can be seen as an approximate slowly rotating solution. In our numerical calculation, we restrict the values of $\ell$ and ${\tilde{a}}$ to [0,  0.2] and [0,  0.3], respectively, to ensure that $\delta \approx 0$.

Appendix B: Coefficients of the scalar perturbation

The coefficients appearing in equation (29) are listed in this Appendix and the relation $\Lambda =l(l+1)$ is also used in this appendix. The coefficients are as follows:

$$\begin{aligned} A_l&=a^2\left( \frac{4+\Lambda +r^2\left( \mu ^2-\omega ^2\right) }{2r^4(r-2M)}-\frac{M\left( 4+\Lambda +r^2\mu ^2\right) }{r^5(r-2M)}\right. \nonumber \\&\quad \left. +\frac{m^2(1+\ell )}{r^4(r-2M)}-\frac{4M^2(1+\ell )\omega ^2}{r^3(r-2M)^2}-\frac{3(r-3M)}{r^6(1+\ell )}\right) \Psi _l\nonumber \\&\quad +\frac{2M\left( \Lambda +r^2\mu ^2\right) }{r^3(r-2M)}\Psi _l-\frac{\Lambda +r^2\left( \mu ^2-\omega ^2\right) }{r^2(r-2M)}\Psi _l\nonumber \\&\quad +\frac{2M}{r^4(1+\ell )}\left( r\Psi '_l-\Psi _l\right) +\frac{r-2M}{r^2(1+\ell )}\Psi ''_l\nonumber \\&\quad +\frac{a^2(r+2M+2r\ell )}{2r^4(1+\ell )}\Psi ''_l-\frac{a^2(5M+2r\ell )}{r^5(1+\ell )}\Psi '_l\nonumber \\&\quad -\frac{4mMa\sqrt{1+\ell }\omega }{r^3(r-2M)}\Psi _l, \end{aligned}$$

(B1)

$$\begin{aligned} D_l&=\frac{2M}{r^6}\Psi _l-\frac{(1+\ell )\left[ 2M\left( \Lambda +r^2\omega ^2\right) -r\Lambda \right] }{r^5(r-2M)}\Psi _l\nonumber \\&\quad -\frac{2M}{r^5}\Psi _l'-\frac{r-2M}{r^4}\Psi _l''. \end{aligned}$$

(B2)

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

$$\begin{aligned} A^{(0)}_l&=\frac{i(r-4M)\omega }{r^2(r-2M)(1+\ell )}\textrm{R}_l-\frac{\Lambda +r^2\mu ^2}{r^3}\textrm{P}_l-\frac{i\omega }{r^2}\textrm{S}_l\nonumber \\&\quad +\frac{i\omega }{r(1+\ell )}\textrm{R}'_l+\frac{r-2M}{r^2(1+\ell )}\textrm{P}''_l+\frac{2mMa\sqrt{1+\ell }}{r^3(r-2M)}\nonumber \\&\quad \times \left[ \frac{i}{r(1+\ell )}\left( \textrm{R}_l-\frac{r-2M}{\Lambda }\textrm{S}'_l \right) -\omega \textrm{P}_l-\frac{ir\omega ^2}{\Lambda }\textrm{S}_l \right] \nonumber \\&\quad +\frac{m^2a^2(1+\ell )}{r^4(r-2M)}\left( \textrm{P}_l+\frac{i\omega r}{\Lambda }\textrm{S}_l\right) +\frac{ia^2\omega }{r^3(1+\ell )}\textrm{R}'_l\nonumber \\&\quad +\frac{2a^2(r+M)}{r^6}\left( \textrm{P}_l-r\textrm{P}'_l+\frac{r^3}{2(r+M)}\textrm{P}''_l\right) \nonumber \\&\quad +\frac{i\omega a^2}{r-2M}\left( \frac{4M^2}{r^4(r-2M)} -\frac{1}{r^3(1+\ell )}\right) \textrm{R}_l, \end{aligned}$$

(C1)

$$\begin{aligned} A^{(1)}_l&=\frac{2M\left( \Lambda +r^2\mu ^2\right) -r\Lambda -r^3\left( \mu ^2-\omega ^2\right) }{r^2(r-2M)^2}\textrm{R}_l+\frac{1}{r^2}\textrm{S}'_l\nonumber \\&\quad +\frac{i \omega }{r(r-2M)}\textrm{P}_l-\frac{i\omega }{r-2M}\textrm{P}'_l-\frac{2iaMm\sqrt{1+\ell }}{r^2(r-2M)\Lambda }\nonumber \\&\quad \times \left[ \frac{\Lambda }{r^2}\textrm{P}_l-\frac{2i\Lambda \omega }{r-2M}\textrm{R}_l-\frac{\Lambda }{r}\textrm{P}'_l+i\omega \textrm{S}'_l \right] \nonumber \\&\quad +a^2\left( \frac{r\omega ^2-\mu ^2(r-2M)}{r^2(r-2M)^2}-\frac{\Lambda }{r^4(r-2M)} \right) \textrm{R}_l\nonumber \\&\quad +\frac{(1+\ell )a^2}{r^4(r-2M)^2}\left[ \left( r\Lambda -\frac{4r^2M^2+r^4}{r-2M} \right) \textrm{R}_l\right. \nonumber \\&\quad +\left. \left( rm^2+r^3\mu ^2\right) \textrm{R}_l+4iM^2\omega \left( r\textrm{P}'_l-\textrm{P}_l\right) \right. \nonumber \\&\quad \left. -\frac{r(r-2M)m^2}{\Lambda }\textrm{S}'_l\right] \end{aligned}$$

(C2)

$$\begin{aligned} A^{(L)}_l&=-\frac{2aMm\sqrt{1+\ell }}{\Lambda (r-2M)r^3}\left( i\Lambda \textrm{P}_l+r\omega \textrm{S}_l \right) +\frac{a^2(1+\ell )}{r^3(r-2M)}\nonumber \\&\quad \times \left( \frac{m^2}{\Lambda }\textrm{S}_l-\frac{4iM^2\omega }{r-2M}\textrm{P}_l\right) +\frac{a^2}{r^4(1+\ell )}\left( r\textrm{R}'_l-\textrm{R}_l\right) \nonumber \\&\quad +\frac{i \omega }{r-2M}\textrm{P}_l+\frac{\textrm{R}_l+r\textrm{R}'_l}{r^2(1+\ell )}-\frac{1}{r^2}\textrm{S}_l, \end{aligned}$$

(C3)

$$\begin{aligned} {\tilde{A}}^{(0)}_l&=\frac{4aM\sqrt{1+\ell }}{r^5}\textrm{Q}_l-\frac{4a^2Mm(1+\ell )\omega }{\Lambda r^5}\textrm{Q}_l, \end{aligned}$$

(C4)

$$\begin{aligned} {\tilde{A}}^{(1)}_l&=-\frac{2ia^2m(1+\ell )}{\Lambda r^4}\textrm{Q}'_l, \end{aligned}$$

(C5)

$$\begin{aligned} {\tilde{A}}^{(L)}_l&=0, \end{aligned}$$

(C6)

$$\begin{aligned} B^{(0)}_l&=\frac{2aM\sqrt{1+\ell }\omega ^2}{r^2(r-2M)\Lambda }\textrm{Q}_l+\frac{2aM}{r^4\sqrt{1+\ell }\Lambda }\textrm{Q}'_l\nonumber \\&\quad -\frac{a^2m(1+\ell )\omega }{r^3(r-2M)\Lambda }\textrm{Q}_l, \end{aligned}$$

(C7)

$$\begin{aligned} B^{(1)}_l&=\frac{2i aM\sqrt{1+\ell }\omega }{\Lambda (r-2M)r^2}\textrm{Q}'_l-\frac{ia^2m(1+\ell )}{\Lambda (r-2M)r^3}\textrm{Q}'_l, \end{aligned}$$

(C8)

$$\begin{aligned} B^{(L)}&= -\frac{2iaM\sqrt{1+\ell }\omega }{\Lambda (r-2M)r^2}\textrm{Q}_l+\frac{ia^2m(1+\ell )}{\Lambda (r-2M)r^3}\textrm{Q}, \end{aligned}$$

(C9)

$$\begin{aligned} {\tilde{B}}^{(0)}_l&=\frac{4a^2M(1+\ell )\left( \Lambda \textrm{P}_l+ir\omega \textrm{S}_l\right) }{\Lambda r^6}, \end{aligned}$$

(C10)

$$\begin{aligned} {\tilde{B}}^{(1)}_l&=\frac{2a^2(1+\ell )\left[ \Lambda \textrm{R}_l-(r-2M)\textrm{S}'_l\right] }{\Lambda (r-2M)r^4}, \end{aligned}$$

(C11)

$$\begin{aligned} {\tilde{B}}^{(L)}_l&=0, \end{aligned}$$

(C12)

$$\begin{aligned} D^{(0)}_l&= \frac{a^2}{(r-2M)r^6}\left[ (2M-r)\left( 6M-\Lambda r(1+\ell ) \right) \textrm{P}_l\right. \nonumber \\&\quad \left. +r\left( i(10M-r)r\omega \textrm{R}_l-(2M-r)\left( ir (1+\ell )\omega \textrm{S}_l\right. \right. \right. \nonumber \\&\quad \left. \left. \left. +6M\textrm{P}'_l+r\left( -ir\omega \textrm{R}'_l+(2M-r)\textrm{P}''_l \right) \right) \right) \right] , \end{aligned}$$

(C13)

$$\begin{aligned} D^{(1)}_l&=\frac{a^2(1+\ell )}{r^4(r-2M)^2}\left[ \left( \Lambda (r-2M)-2Mr^2\omega ^2 \right) \textrm{R}_l\right. \nonumber \\&\quad \times (r-2M)\left( 2iMr\omega \textrm{P}'_l-(r-2M)\textrm{S}'_l \right) \nonumber \\&\quad \left. +2iM(2M-r)\omega \textrm{P}_l \right] , \end{aligned}$$

(C14)

$$\begin{aligned} D^{(L)}_l&=-\frac{a^2}{(r-2M)r^4}\left[ 2iMr(1+\ell )\omega \textrm{P}_l\right. \nonumber \\&\quad \left. +(r-2M)\left( \textrm{R}_l-(1+\ell )\textrm{S}_l+r\textrm{R}'_l \right) \right] , \end{aligned}$$

(C15)

$$\begin{aligned} \alpha _l&=\frac{(1+\ell )^{-1}}{(2M-r)r^2\Lambda }\big [ir^2(1+\ell )\Lambda \omega \textrm{P}_l-(r-2M)\Lambda \textrm{R}_l\nonumber \\&\quad +r^2(1+\ell )\left( r\mu ^2-r\omega ^2-2M\mu ^2 \right) \textrm{S}_l+r\Lambda (r-2M)\textrm{R}'_l \nonumber \\&\quad -(r-2M)\left( 2M\textrm{S}'_l+r(r-2M)\textrm{S}''_l \right) \big ] \nonumber \\&\quad +\frac{2aMm\sqrt{1+\ell }}{\Lambda (r-2M)r^3}\left( i\Lambda \textrm{P}_l-r\omega \textrm{S}_l \right) +\frac{a^2(1+\ell )^{-1}}{r^4(r-2M)^2\Lambda }\nonumber \\&\quad \times \left[ 2i\Lambda Mr^2(1+\ell )^2\omega \textrm{P}_l-3\Lambda (r-2M)^2\ell \textrm{R}_l\right. \nonumber \\&\quad \left. +(r-2M)^2\left( \Lambda r\ell \textrm{R}'_l-2M(1+\ell )\left( 3\textrm{S}'_l-r\textrm{S}''_l \right) \right) \right. \nonumber \\&\quad \left. -2Mr^3(1+\ell )^2\omega ^2\textrm{S}_l \right] , \end{aligned}$$

(C16)

$$\begin{aligned} \beta _l&=\left( \frac{1}{r^2}+\frac{\mu ^2}{\Lambda }-\frac{r\omega ^2}{\Lambda (r-2M)} \right) \textrm{Q}_l-\frac{2M}{r^2\Lambda (1+\ell )}\textrm{Q}'_l\nonumber \\&\quad -\frac{(r-2M)}{r\Lambda (1+\ell )}\textrm{Q}''_l+\frac{2aMm\sqrt{1+\ell }\omega }{\Lambda (r-2M)r^2}\textrm{Q}_l\nonumber \\&\quad -\frac{(1+\ell )a^2}{r^4}\textrm{Q}_l+\frac{2Ma^2(1+\ell )\omega ^2}{r(r-2M)^2\Lambda }\textrm{Q}_l\nonumber \\&\quad +\frac{2Ma^2}{r^3\Lambda }\left( \frac{3}{r}\textrm{Q}'_l-\textrm{Q}''_l\right) , \end{aligned}$$

(C17)

$$\begin{aligned} \zeta _l&=\frac{6aM}{r^4\sqrt{1+\ell }}\left( \frac{ir^2\omega }{r-2M}\textrm{R}_l-\textrm{P}_l+r\textrm{P}'_l \right) +\frac{a^2m}{r^4\Lambda }\nonumber \\&\quad \times \left[ \frac{2Mr(1+\ell )\Lambda \omega }{r-2M} \textrm{P}_l-\frac{i(8M-r)\Lambda }{r-2M}\textrm{R}_l-ir\Lambda \textrm{R}'_l\right. \nonumber \\&\quad \left. +\frac{2iMr^2(1+\ell )\omega ^2}{r-2M}\textrm{S}_l+i\left( 8M\textrm{S}'_l+(r-2M)r\textrm{S}''_l\right) \right] , \end{aligned}$$

(C18)

$$\begin{aligned} \eta _l&=\frac{2iaM\sqrt{1+\ell }\omega }{(2M-r)r^2}\textrm{Q}_l+\frac{ia^2m}{r^4\Lambda }\bigg [\frac{2M(1+\ell )\left( \Lambda +r^2\omega ^2\right) }{r-2M}\textrm{Q}_l\nonumber \\&\quad -2\left( r-3M\right) \textrm{Q}'_l+r(r-2M)\textrm{Q}''_l\bigg ], \end{aligned}$$

(C19)

$$\begin{aligned} \rho _l&=\frac{a^2}{r^4(r-2M)\Lambda }\left[ 2Mr(1+\ell )\omega \left( r\omega \textrm{S}_l-i\Lambda \textrm{P}_l\right) \right. \nonumber \\&\quad \left. -2(r-3M)(r-2M)\textrm{S}'_l+r(r-2M)^2\textrm{S}''_l\right. \nonumber \\&\quad \left. +\Lambda (r-2M)\left( 3\textrm{R}_l-r\textrm{R}'_l\right) \right] , \end{aligned}$$

(C20)

$$\begin{aligned} \lambda _l&=-\frac{4a^2M(1+\ell )}{r^5}\textrm{Q}_l, \end{aligned}$$

(C21)

$$\begin{aligned} \sigma _l&=0, \end{aligned}$$

(C22)

$$\begin{aligned} \gamma _l&=\frac{a^2}{r^4(r-2M)\Lambda }\left[ (1+\ell )\left( r\Lambda -2M(\Lambda +r^2\omega ^2)\right) \textrm{Q}_l\right. \nonumber \\&\quad \left. -8M(r-2M)\textrm{Q}'_l-r(r-2M)^2\textrm{Q}''_l \right] . \end{aligned}$$

(C23)

After the angular parts are separated, Eqs. (45)–(47) can be rewritten as the pure radial equations, which are given by

$$\begin{aligned}&A^{(I)}_l+{\mathcal {Q}}^2_{l+1}\left[ D^{(I)}_l+l{\tilde{B}}^{(I)}_l \right] +{\mathcal {Q}}^2_l\left[ D^{(I)}_l-(l+1){\tilde{B}}^{(I)}_l\right] \nonumber \\&\quad +{\mathcal {Q}}_l\left[ {\tilde{A}}^{(I)}_{l-1}+(l-1)B^{I}_{l-1} \right] +{\mathcal {Q}}_{l+1}\left[ {\tilde{A}}^{(I)}_{l+1}-(l+2)B^{(I)}_{l+1} \right] \nonumber \\&\quad +{\mathcal {Q}}_{l-1}{\mathcal {Q}}_l\left[ D^{(I)}_{l-2}+(l-2){\tilde{B}}^{(I)}_{l-2} \right] \nonumber \\&\quad +{\mathcal {Q}}_{l+2}{\mathcal {Q}}_{l+1}\left[ D^{(I)}_{l+2}-(l+3){\tilde{B}}^{(I)}_{l+2} \right] =0, \end{aligned}$$

(C24)

$$\begin{aligned}&{\mathcal {Q}}_{l+1}\left[ l \eta _{l+1}+i m\left( (l+2)\gamma _{l+1}-\lambda _{l+1}\right) \right] \nonumber \\&\quad -{\mathcal {Q}}_{l-1}{\mathcal {Q}}_l(l+1)\left[ (l-2)\rho _{l-2}+\sigma _{l-2} \right] \nonumber \\&\quad +{\mathcal {Q}}_{l+2}{\mathcal {Q}}_{l+1}l\left[ -(l+3)\rho _{l+2}+\sigma _{l+2} \right] \nonumber \\&\quad +\Lambda \alpha _l-i m\zeta _l+{\mathcal {Q}}^2_{l+1}l\left[ l\rho _l+\sigma _l \right] \nonumber \\&\quad +{\mathcal {Q}}^2_l(l+1)\left[ (l+1)\rho _l-\sigma _l\right] \nonumber \\&\quad +{\mathcal {Q}}_l\left[ -(l+1)\eta _{l-1}-i m\left( (l-1)\gamma _{l-1}+\lambda _{l-1}\right) \right] =0, \end{aligned}$$

(C25)

$$\begin{aligned}&{\mathcal {Q}}_{l+1}\left[ l \zeta _{l+1}-i m\left( (l+2)\rho _{l+1}-\sigma _{l+1}\right) \right] \nonumber \\&\quad -{\mathcal {Q}}_{l-1}{\mathcal {Q}}_l(l+1)\left[ (l-2)\gamma _{l-2}+\lambda _{l-2} \right] \nonumber \\&\quad +{\mathcal {Q}}_{l+2}{\mathcal {Q}}_{l+1}l\left[ -(l+3)\gamma _{l+2}+\lambda _{l+2} \right] \nonumber \\&\quad +\Lambda \beta _l+i m\eta _l+{\mathcal {Q}}^2_{l+1}l\left[ l\gamma _l+\lambda _l \right] \nonumber \\&\quad +{\mathcal {Q}}^2_l(l+1)\left[ (l+1)\gamma _l-\lambda _l\right] \nonumber \\&\quad +{\mathcal {Q}}_l\left[ -(l+1)\zeta _{l-1}+i m\left( (l-1)\rho _{l-1}+\sigma _{l-1}\right) \right] =0. \end{aligned}$$

(C26)

The source term corresponding to the coupled system described by Eqs. (56), (57) is

$$\begin{aligned} {\mathcal {S}}^{\text {po}}_1&=(1+\ell )^{\frac{3}{2}}\frac{6 i a M{\mathcal {F}}\omega }{\Lambda r^3}\left[ l{\mathcal {Q}}_{l+1}\textrm{Q}_{l+1}-(l+1){\mathcal {Q}}_l\textrm{Q}_{l-1} \right] \nonumber \\&\quad +\sqrt{1+\ell }\frac{2 a M m}{\Lambda r^5\omega }\left[ \Lambda \left( 2(1+\ell )r^2\omega ^2+3 {\mathcal {F}}^2\right) \textrm{R}_l\right. \nonumber \\&\quad \left. +3{\mathcal {F}}\left( r \Lambda {\mathcal {F}}~(\textrm{R}_l)'-(1+\ell )\left( \Lambda {\mathcal {F}}+r^2\omega ^2\right) \textrm{S}_l \right) \right] , \end{aligned}$$

(C27)

$$\begin{aligned} {\mathcal {S}}^{\text {po}}_2&=\sqrt{1+\ell }\frac{2a Mm}{r^5\omega }\left[ 2(1+\ell )r^2\omega ^2\textrm{S}_l+3r{\mathcal {F}}^2\textrm{S}_l'\right. \nonumber \\&\quad \left. -3{\mathcal {F}}\left( \Lambda +r^2\mu ^2 \right) \textrm{R}_l \right] , \end{aligned}$$

(C28)

and the corresponding source term of Eq. (58) is

$$\begin{aligned} {\mathcal {S}}^{\text {ax}}= -\frac{6i aM{\mathcal {F}}}{\sqrt{1+\ell }r^5\omega }\left[ (l+1){\mathcal {Q}}_l\psi _{l-1}-l{\mathcal {Q}}_{l+1}\psi _{l+1} \right] , \end{aligned}$$

(C29)

where we have defined the polar function

$$\begin{aligned} \psi _l=(1+\ell )\left[ (\Lambda +r^2\mu ^2)\textrm{R}_l-(r-2M)\textrm{S}_l'\right] . \end{aligned}$$

(C30)

Note that this term is similar to Eq. (25) in Ref. [28]. As expected, the axial perturbation $\textrm{Q}_l$ is coupled to the polar functions with $l\pm 1$.

Appendix D: Recurrence coefficients of gravitational perturbation

In this appendix, we present the explicit expressions of the recurrence coefficients when we used the continued fraction method. The recurrence coefficients in Eq. (87) are

$$\begin{aligned} n_1{\mathcal {C}}_{1,0}&=\tilde{\ell }m{\tilde{a}}[3i+\Lambda \omega (1-3i\omega )]\nonumber \\&\quad -i\Lambda \omega [\Lambda -3-4\omega (i+2\omega )], \end{aligned}$$

(D1)

$$\begin{aligned} n_2{\mathcal {C}}_{2,0}&=\tilde{\ell }m{\tilde{a}}[30i+\Lambda (1-2i\omega )\omega ]\nonumber \\&\quad \times i\Lambda \omega [3+4\omega (i+\omega )], \end{aligned}$$

(D2)

$$\begin{aligned} n_2{\mathcal {C}}_{2,1}&=3i\tilde{\ell }m{\tilde{a}}[\Lambda \omega (i+\omega )-1]\nonumber \\&\quad +i\Lambda \omega [1+\Lambda -4\omega (3i+2\omega )], \end{aligned}$$

(D3)

$$\begin{aligned} n_3{\mathcal {C}}_{3,0}&=72i\tilde{\ell }m{\tilde{a}}, \end{aligned}$$

(D4)

$$\begin{aligned} n_3{\mathcal {C}}_{3,1}&=4\Lambda (2-i\omega )\omega ^2+2i\tilde{\ell }m{\tilde{a}}[\Lambda \omega (i+\omega )-15], \end{aligned}$$

(D5)

$$\begin{aligned} n_3{\mathcal {C}}_{3,2}&=\tilde{\ell }m{\tilde{a}}[3i+\Lambda (5-3i\omega )\omega ]\nonumber \\&\quad -i\Lambda \omega [9+\Lambda -4\omega (5i+2\omega )], \end{aligned}$$

(D6)

$$\begin{aligned} n_4{\mathcal {C}}_{4,0}&=66i\tilde{\ell }m{\tilde{a}}, \end{aligned}$$

(D7)

$$\begin{aligned} n_4{\mathcal {C}}_{4,1}&=-72i\tilde{\ell }m{\tilde{a}}, \end{aligned}$$

(D8)

$$\begin{aligned} n_4{\mathcal {C}}_{4,2}&=\tilde{\ell }m{\tilde{a}}[30i+\Lambda (3-2i\omega )\omega ]\nonumber \\&\quad +i\Lambda \omega [4\omega (3i+\omega )-5], \end{aligned}$$

(D9)

$$\begin{aligned} n_4{\mathcal {C}}_{4,3}&=i\Lambda \omega [21+\Lambda -4\omega (7i+2\omega )]\nonumber \\&\quad +i\tilde{\ell }m{\tilde{a}}[\Lambda \omega (7i+3\omega )-3], \end{aligned}$$

(D10)

where

$$\begin{aligned} n_1&=\Lambda \omega (\tilde{\ell }m{\tilde{a}}-2\omega -i), \end{aligned}$$

(D11)

$$\begin{aligned} n_2&=2\Lambda \omega (2i-\tilde{\ell }m{\tilde{a}}+2\omega ), \end{aligned}$$

(D12)

$$\begin{aligned} n_3&= 3\Lambda \omega (\tilde{\ell }m{\tilde{a}}-2\omega -3i), \end{aligned}$$

(D13)

$$\begin{aligned} n_4&=4\Lambda \omega (4i-\tilde{\ell }m{\tilde{a}}+2\omega ). \end{aligned}$$

(D14)

The recurrence coefficients in Eq. (88) are

$$\begin{aligned} \alpha _n&=4\Lambda (1+n)(1+n+i\tilde{\ell }m{\tilde{a}}-2i\omega )\omega , \end{aligned}$$

(D15)

$$\begin{aligned} \beta _m&= 4\Lambda \omega [3-\Lambda -2n^2-n(2-8i\omega )+4i\omega +8\omega ^2]\nonumber \\&\quad +4\tilde{\ell }m{\tilde{a}}[3-\Lambda \omega (i+2in+3\omega )], \end{aligned}$$

(D16)

$$\begin{aligned} \gamma _m&=4\Lambda \omega (n^2-4-4in\omega -4\omega ^2)\nonumber \\&\quad +4\tilde{\ell }m{\tilde{a}}[i\Lambda n \omega +2(\Lambda \omega ^2-15)], \end{aligned}$$

(D17)

$$\begin{aligned} \sigma _n&=288\tilde{\ell }m{\tilde{a}}, \end{aligned}$$

(D18)

$$\begin{aligned} \tau _n&=-264\tilde{\ell }m{\tilde{a}}, \end{aligned}$$

(D19)

$$\begin{aligned} \delta _n&=84\tilde{\ell }m{\tilde{a}}. \end{aligned}$$

(D20)

Appendix E: QNM frequency tables

See Tables 1, 2 and 3.

Table 1 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 2 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 3 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