# Universal horizons and black hole spectroscopy in gravitational theories
with broken Lorentz symmetry

Chao Zhanga,b,c,d<EMAIL_ADDRESS><EMAIL_ADDRESS>Anzhong
Wange111Corresponding author<EMAIL_ADDRESS>Corresponding author
Tao Zhua,b<EMAIL_ADDRESS>a Institute for theoretical physics and
cosmology, Zhejiang University of Technology, Hangzhou, 310032, China
b United Center for Gravitational Wave Physics (UCGWP), Zhejiang University of
Technology, Hangzhou, 310032, China
c College of Information Engineering, Zhejiang University of Technology,
Hangzhou, 310032, China
d Department of Physics, Faculty of Science, Tokyo University of Science, 1-3,
Kagurazaka, Shinjuku-ku, Tokyo 162-8601, Japan
eGCAP-CASPER, Physics Department, Baylor University, Waco, Texas, 76798-7316,
USA

###### Abstract

The violation of Lorentz invariance (LI) in gravitational theories, which
allows superluminal propagations, dramatically alters the causal structure of
the spacetime and modifies the notion of black holes. In particular, instead
of metric horizons, now universal horizons define the boundaries of black
holes, within which a particle cannot escape to spatial infinities even with
an infinitely large speed. Then, a natural question is how the quasi-normal
modes (QNMs) of a black hole are modified, if one considers the universal
horizon as its causal boundary. In this paper, we study in detail this problem
in Einstein-aether theory, a vector-tensor theory that violates LI but yet is
self-consistent and satisfies all observations carried out so far.
Technically, this poses several challenges, including singularities of the
perturbation equations across metric horizons and proper identifications of
ingoing modes at universal horizons. After overcoming these difficulties, we
show that the QNMs of the Schwarzschild black hole, also a solution of
Einstein-aether theory, consist of two parts, the metric and aether parts. The
QNMs of the metric perturbations are quite similar to those obtained in
general relativity and are consistent with current observations of
gravitational waves. But the ones from aether perturbations are different, and
our numerical studies indicate that they are even not stable. The latter is
consistent with our previous studies, which showed that the stealth
Schwarzschild black hole suffers a Laplacian instability along the angular
direction. The method and techniques developed in this paper can be applied to
the studies of QNMs in other theories of gravity with broken LI.

##  I Introduction

The detection of the first gravitational wave (GW) from the coalescence of two
massive black holes (BHs) by advanced LIGO marked the beginning of a new era —
the GW astronomy Ref1 . Following this observation, about 90 GW events have
been identified by the LIGO/Virgo/KAGRA (LVK) scientific collaborations (see,
e.g., GWs ; GWs19a ; GWs19b ; GWsO3b ). In the future, more ground- and space-
based GW detectors will be constructed Moore2015 ; Gong:2021gvw , which will
enable us to probe signals with a much wider frequency band and larger
distances. This triggered the interests on the quasi-normal mode (QNM) of
black holes, as GWs emitted in the ringdown phase can be considered as the
linear combination of these individual modes Berti2009 ; Berti18 .

From the classical point of view, QNMs of black holes are eigenmodes of
dissipative systems. The information contained in QNMs provide the keys in
revealing whether BHs are ubiquitous in our universe, and more important
whether general relativity (GR) is the correct theory to describe gravity even
in the strong field regime Berti18b . Basically, a QNM frequency $\omega$
contains two parts, the real and imaginary parts. Its real part gives the
frequency of vibration while its imaginary part provides the damping time.

In GR, according to the no-hair theorem, an isolated and stationary BH is
completely characterized by only three quantities, mass, angular momentum and
electric charge. Astrophysically, we expect BHs to be neutral, so they are
uniquely described by the Kerr solution. Then, the QNM frequencies and damping
times will depend only on the mass and angular momentum of the finally formed
BH. Clearly, to extract physics from the ringdown phase, at least two QNMs are
needed. This will require the signal-to-noise ratio (SNR) to be of the order
100 Berti18 . Although such high SNRs are not achievable right now, it has
been shown that they may be achievable once the advanced LIGO, Virgo and KAGRA
reach their fully designed sensitivities. In any case, it is certain that they
will be detected by the ground-based third-generation detectors, such as
Cosmic Explorer CE  and the Einstein Telescope ET , as well as the space-based
detectors, including LISA LISA , TianQin Liu2020 , Taiji Taiji2 , and DECIGO
DECIGO .

QNMs in GR have been studied extensively Chandra92 , including scalar, vector
and tensor perturbations Iyer1987 . Such calculations have been extended from
the Schwarzschild BH to other more general cases, e.g., the Kerr BH Det1980 ;
Seidel1990 . In this procedure, several different techniques of computations
of QNMs were developed. For instance, the Wentzel-Kramers-Brillouin (WKB)
approach Will1985 ; Will1987 ; Konoplya2003 ; Jerzy2017 , the finite
difference method (FDM) XinLi2020 , the continued fraction method Leaver1985 ,
the shooting method Chandra1975 ; Doneva2010 , the matrix method Kai2017 , and
so on Kono2011 ; Gund1994 ; Bin2004 . Some of these methods have been also
applied to modified theories of gravity Xin2020 ; Oliver2019 . Additionally,
some special approximations, e.g., the eikonal limit, have also been
extensively explored, see, for example, Huan2012  and references therein.

In this paper, we shall focus on the QNMs of black holes in Einstein-aether
theory ($\ae$-theory) JM01 . Such studies are well motivated. In particular,
the theory is self-consistent, such as free of ghosts and instability Jacobson
, and satisfies all the experimental tests carried out so far OMW18 . Its
Cauchy problem is also well posed SBP19 , and energy is always positive (as
far as the hypersurface-orthogonal aether field is concerned) GJ11 . In
addition, there exist both gravitational waves Foster06 ; Foster07 ;
Yagi:2013qpa ; Yagi14 ; HYY15 ; GHLP18 ; Kai19 ; Zhao19 ; Chao2020 ; OMW21 ;
GHBBCYY  and BHs Eling2006-1 ; Eling2006-2 ; Tamaki2008 ; BS11 ; BJS11 ; Per12
; Wang13 ; Dingq15 ; Chao2020b ; AFJW21 ; CSS21 . It was also shown that
universal horizons can be formed from gravitational collapse of realistic
matter BMSW18 , in addition to the formation of the spin-0 and metric horizons
Garfinkle2007 .

In comparison with other modified theories of gravity CFPS12 , including
scalar-tensor theories and their high-order corrections DL19 , $\ae$-theory
has at least two distinguishable features:

  * •

It is a particular vector-tensor theory in which the vector field is always
timelike. As a result, it always defines a preferred frame and whereby
violates locally the Lorentz invariance (LI). Despite the facts that LI is the
cornerstone of modern physics, and all the experiments carried out so far are
consistent with it Bourgoin21 ; Bars19 ; Shao19 ; Bourgoin17 ; Flowers17 ;
KR11 , violations of LI have been well motivated and extensively studied in
the past several decades, especially from the point of view of quantum gravity
Collin04 ; Mattingly05 ; Liberati13 ; PT14 ; Wang17 .

  * •

BHs in $\ae$-theory are defined in terms of universal horizons (UHs) BS11 ;
Bhattacharyya2016 ; Bhattacharyya2016b , instead of metric horizons (MHs) as
that in GR. This is due to the fact that in $\ae$-theory there exist three
different species of gravitons, spin-0, spin-1, and spin-2, and each of them
move in different speeds JM04 . To avoid the vacuum gravi-Čerenkov radiation,
such as cosmic rays, each of these three species must move with a speed that
is at least no less than the speed of light EMS05 . As a matter of fact,
depending on the choice of the free coupling constants of the theory, they can
be arbitrarily large, and so far no upper limits of these speeds are known
KR11 .

Clearly, now MHs are no longer one-way membranes to these particles, and they
can cross MHs to escape to infinities even initially they are trapped inside
them. On the other hand, a UH is defined as the causal boundary of particles
with arbitrarily large speeds BS11 , and once they are inside it, they can
never cross it to escape to spatial infinities (For a recent review of UHs,
see, e.g., Wang17  and references therein).

With the above remarkable features of $\ae$-theory, it would be very
interesting and important to find new predictions of the theory for the BH
spectroscopy Berti18 , the QNMs mentioned above, which has been extensively
studied in the last couple of years in terms of GWs emitted in the ringdown
phase of binary BHs (BBHs) (for example, see BLT21 ; LCV20 ; FBPF20 ; BFPF20 ;
GST19 ; IGFST19 ; CDV19 ; CN20 ; OC20 ; Shaik20 ; Uchi20 ; BCCK20 ; Cabero20 ;
Maselli20 ; Islam20  and references therein), and found that they are all
consistent with GR within the error bars allowed by the observations of the 90
GW events LVK2022 .

As the first step to study QNMs in $\ae$-theory, recently some of the current
authors studied spherically symmetric static spacetimes and obtained various
numerical solutions Chao2020b , with the choice of the four free parameters
involved in the theory that satisfy all the observational constraints OMW18
222It should be noted that all the solutions obtained previously do not
satisfy these conditions, which were the main motivations for us to revisited
this problem and obtained new (numerical) solutions in Chao2020b .. We also
derived a set of analytic solutions in which the four-parameters satisfy all
the self-consistent and observational constraints obtained in OMW18 , and
showed that they are nothing but the (anti- or) de Sitter Schwarzschild black
holes, in which the aether field are non-trivial but does not contribute to
the curvatures of the spacetime. When the cosmological constant vanishes, the
universal horizon locates precisely at $r=3r_{s}/4$, where $r_{s}$ denotes the
location of the MH.

Lately, we extended our above studies to the stability of such obtained new
black holes against odd parity perturbations, and found that they are stable
provided that TZZW21 ,

$c_{4}=0,$ (1.1)

in addition to the constraint OMW18 ,

$c_{13}\equiv c_{1}+c_{3}=0,$ (1.2)

obtained from the observations of GW170817 GW170817  together with its gamma-
ray burst GRB 170817A GRB170817A , where $c_{i}$’s are the four dimensionless
coupling constants of the theory.

In this paper, we continuously study the odd-parity perturbations, but focus
ourselves on the QNMs of the newly found black holes Chao2020b , which, as
mentioned above, satisfy all the constraints of the theory OMW18 . In doing
so, several technical questions raise. One of them is connected with the inner
boundary conditions, as now the inner boundary is the UH, which is always
inside the MH, at which the timelike Killing vector becomes spacelike. Then,
the notion of “in”- and “out”-going waves must be define properly, before
imposing the condition that there are no out-going waves at the UH. In
addition, the field equations of perturbations are usually singular at the MH,
and proper smoothness conditions must be imposed. In GR, the perturbation
equations are also singular at MHs, but such smoothness conditions are nicely
avoided, as now the MH is the inner boundary, and the smoothness conditions
are simply replaced by the inner boundary conditions.

To overcome the above issues, we first choose the ingoing Eddington-
Finkelstein (EF) coordinates, so the background is smoothly across the MH, and
more important, the ingoing and outgoing waves are well defined over the whole
spacetime, including the region between the UH and MH. Then, we carefully
study the structure of the field equations near the MH, and find the
smoothness conditions that we must impose in order to assure that the field
equations hold across the MH. Once the above issues are clarified, we first
impose the purely ingoing wave conditions at the UH and purely outgoing wave
conditions at the spatial infinity, and then solve the field equations of the
odd-parity perturbations by shooting method, respectively, from the UH and the
spatial infinity toward the MH. Then, at the MH we impose the smoothness
conditions, which are satisfied only for particular choices of the QNM
frequency $\omega$, whereby we read off the QNMs of the odd-parity
perturbations.

To show the above explicitly, in this paper we shall take the Schwarzschild
black hole solution as the background, since it is also a solution of the
$\ae$-theory Chao2020b . The main reasons in doing are two-folds: i) First, in
such a background, the mathematics involved will be much simple, from which we
can construct clearly the boundary conditions at UHs and at spatial
infinities, as well as the smoothness conditions across MHs. ii) Second, we
can understand deeply the nature of the instability of this solution found in
TZZW21  from the point of view of QNMs, as in TZZW21  only its stability
against large angular perturbations was studied.

Specifically, the paper is organized as follows. In Section II we first
provide a brief introduction to $\ae$-theory, and then give various
definitions of horizons, including the particle, metric and universal
horizons, which are important to the studies of QNMs of black holes, because
the fact that superluminal motions exist in the theory, and the inner
boundaries of black holes now are the universal horizons. In Section III, we
consider the odd-parity perturbations, and construct the master equations,
respectively, for the metric and aether field perturbations. In Section IV we
first work out step by step the boundary conditions for the master equation of
the metric perturbations at the universal horizon as well as at the spatial
infinity. Then, we find out the smoothness conditions across the metric
horizon. Afterwards, combining the Chandrasekhar-Detweiler method Chandra1975
with the shooting method (one side is from the UH to the MH, while the other
side is from the spatial infinity to the MH [cf. Section IV.D]), we solve the
master equation and find that only for particular choices of the mode
$\omega$, can the smoothness conditions at the MH be satisfied, whereby
various values of $\omega$ are found and given explicitly in Table 2 for
$l=2,3,4,10$, respectively. We also compare such a spectrum with the
corresponding one given in GR, and find that the differences between them
cannot be distinguished by current observations of GWs, although it is well
within the detectability of the next generation detectors. In Section V,
closely following what is done for the metric perturbations in Section IV, we
study the aether perturbations and find various modes for the aether field
with $l=2,10,100$, given, respectively, in Tables 4 and 5, in which some
unstable modes are identified. This strongly indicates the instability of the
aether field, and is also consistent with what was found in TZZW21 . However,
due to our numerical errors of the current paper, the existence of such
unstable modes is not exclusive. In particular, in the metric perturbations
presented in Section IV, no such instability is found. The paper is ended with
Section VI, in which we derive our main conclusions and present some remarks
on the future investigations.

Through out the paper, we shall adopt the unit system so that $c=G_{N}=1$,
where $c$ is the speed of light while $G_{N}$ stands for the Newtonian
constant. All the Greek letters run from 0 to 3.

##  II Einstein-aether theory

In Einstein-aether theory ($\ae$-theory), the fundamental variables of the
gravitational sector are Jacobson ,

$\left(g_{\mu\nu},u^{\mu},\lambda\right),$ (2.1)

where $g_{\mu\nu}$ is the four-dimension metric of the spacetime with the
signature $(-,+,+,+)$, $u^{\mu}$ is the aether field, and $\lambda$ is a
Lagrangian multiplier, which guarantees that the aether field is always
timelike and unity. Then, the general action of the theory is given by,

$S=S_{\ae}+S_{m},$ (2.2)

where $S_{m}$ denotes the action of matter, and $S_{\ae}$ the gravitational
action of the $\ae$-theory, given, respectively, by

$\displaystyle S_{\ae}$ $\displaystyle=$ $\displaystyle\frac{1}{16\pi
G_{\ae}}\int{\sqrt{-g}\;d^{4}x\Big{[}{\cal{L}}_{\ae}\left(g_{\mu\nu},u^{\alpha},c_{i}\right)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+{\cal{L}}_{\lambda}\left(g_{\mu\nu},u^{\alpha},\lambda\right)\Big{]},$
$\displaystyle S_{m}$ $\displaystyle=$
$\displaystyle\int{\sqrt{-g}\;d^{4}x\Big{[}{\cal{L}}_{m}\left(g_{\mu\nu},u^{\alpha};\hat{\psi}\right)\Big{]}}.$
(2.3)

Here $\hat{\psi}$ collectively denotes the matter fields, and $g$ is the
determinant of $g_{\mu\nu}$, and

$\displaystyle{\cal{L}}_{\lambda}$ $\displaystyle\equiv$
$\displaystyle\lambda\left(g_{\alpha\beta}u^{\alpha}u^{\beta}+1\right),$
$\displaystyle{\cal{L}}_{\ae}$ $\displaystyle\equiv$ $\displaystyle
R(g_{\mu\nu})-M^{\alpha\beta}_{~{}~{}~{}~{}\mu\nu}\left(D_{\alpha}u^{\mu}\right)\left(D_{\beta}u^{\nu}\right),$
(2.4)

where $D_{\mu}$ denotes the covariant derivative with respect to $g_{\mu\nu}$,
$R$ is the Ricci scalar, and $M^{\alpha\beta}_{~{}~{}~{}~{}\mu\nu}$ is defined
as

$\displaystyle M^{\alpha\beta}_{~{}~{}~{}~{}\mu\nu}\equiv
c_{1}g^{\alpha\beta}g_{\mu\nu}+c_{2}\delta^{\alpha}_{\mu}\delta^{\beta}_{\nu}+c_{3}\delta^{\alpha}_{\nu}\delta^{\beta}_{\mu}-c_{4}u^{\alpha}u^{\beta}g_{\mu\nu},$

with $\delta_{\mu\nu}$ representing the Kronecker delta. Note that here we
assume that matter fields couple not only to $g_{\mu\nu}$ but also to the
aether field $u^{\mu}$. However, in order to satisfy the severe observational
constraints, such a coupling in general is assumed to be absent Jacobson .

The four coupling constants $c_{i}$’s are all dimensionless, and $G_{\ae}$ is
related to the Newtonian constant $G_{N}$ via the relation CL04 ,

$G_{N}=\frac{G_{\ae}}{1-\frac{1}{2}c_{14}},$ (2.6)

where $c_{ij}\equiv c_{i}+c_{j}$.

The variations of the total action, respectively, with respect to
$g_{\mu\nu}$, $u^{\mu}$ and $\lambda$ yield the field equations Chao2020b ,

$\displaystyle R^{\mu\nu}-\frac{1}{2}g^{\mu\nu}R-S^{\mu\nu}$ $\displaystyle=$
$\displaystyle 8\pi G_{\ae}T^{\mu\nu},$ (2.7) $\displaystyle\AE_{\mu}$
$\displaystyle=$ $\displaystyle 8\pi G_{\ae}T_{\mu},$ (2.8) $\displaystyle
g_{\alpha\beta}u^{\alpha}u^{\beta}$ $\displaystyle=$ $\displaystyle-1,$ (2.9)

where $R^{\mu\nu}$ denotes the Ricci tensor, and

$\displaystyle S_{\alpha\beta}$ $\displaystyle\equiv$ $\displaystyle
D_{\mu}\Big{[}J^{\mu}_{\;\;\;(\alpha}u_{\beta)}+J_{(\alpha\beta)}u^{\mu}-u_{(\beta}J_{\alpha)}^{\;\;\;\mu}\Big{]}$
$\displaystyle+c_{1}\Big{[}\left(D_{\alpha}u_{\mu}\right)\left(D_{\beta}u^{\mu}\right)-\left(D_{\mu}u_{\alpha}\right)\left(D^{\mu}u_{\beta}\right)\Big{]}$
$\displaystyle+c_{4}a_{\alpha}a_{\beta}+\lambda
u_{\alpha}u_{\beta}-\frac{1}{2}g_{\alpha\beta}J^{\delta}_{\;\;\sigma}D_{\delta}u^{\sigma},$
$\displaystyle\AE_{\mu}$ $\displaystyle\equiv$ $\displaystyle
D_{\alpha}J^{\alpha}_{\;\;\;\mu}+c_{4}a_{\alpha}D_{\mu}u^{\alpha}+\lambda
u_{\mu},$ $\displaystyle T^{\mu\nu}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{\sqrt{-g}}\frac{\delta\left(\sqrt{-g}{\cal{L}}_{m}\right)}{\delta
g_{\mu\nu}},$ $\displaystyle T_{\mu}$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{\sqrt{-g}}\frac{\delta\left(\sqrt{-g}{\cal{L}}_{m}\right)}{\delta
u^{\mu}},$ (2.10)

with

$J^{\alpha}_{\;\;\;\mu}\equiv
M^{\alpha\beta}_{~{}~{}~{}~{}\mu\nu}D_{\beta}u^{\nu}\,,\quad a^{\mu}\equiv
u^{\alpha}D_{\alpha}u^{\mu}.$ (2.11)

From Eq.(2.8), we find that

$\lambda=u_{\beta}D_{\alpha}J^{\alpha\beta}+c_{4}a^{2}-8\pi
G_{\ae}T_{\alpha}u^{\alpha},$ (2.12)

where $a^{2}\equiv a_{\lambda}a^{\lambda}$.

It is easy to show that the Minkowski spacetime is a solution of $\ae$-theory,
in which the aether is aligned along the time direction,
${u}_{\mu}=\delta^{0}_{\mu}$. Then, the linear perturbations around the
Minkowski background show that the theory in general possess three types of
excitations, scalar (spin-0), vector (spin-1) and tensor (spin-2) modes, with
their squared speeds given by JM04

$\displaystyle c_{S}^{2}$ $\displaystyle=$
$\displaystyle\frac{c_{123}(2-c_{14})}{c_{14}(1-c_{13})(2+c_{13}+3c_{2})}\,,$
$\displaystyle c_{V}^{2}$ $\displaystyle=$
$\displaystyle\frac{2c_{1}-c_{13}(2c_{1}-c_{13})}{2c_{14}(1-c_{13})}\,,$
$\displaystyle c_{T}^{2}$ $\displaystyle=$ $\displaystyle\frac{1}{1-c_{13}},$
(2.13)

respectively. Here $c_{ijk}\equiv c_{i}+c_{j}+c_{k}$.

Requiring that the theory: (a) be self-consistent, such as free of ghosts; and
(b) satisfies all the observational constraints obtained so far, it was found
that the parameter space of the theory is considerably restricted. In
particular, $c_{14}$, $c_{2}$ and $c_{13}$ are restricted to OMW18  333The
recent studies of the neutron binary systems showed that the PPN parameter
$\alpha_{1}$ is further restricted to $|\alpha_{1}|<10^{-5}$ GHBBCYY , which
is an order of magnitude stronger than the bounds from lunar laser ranging
experiments MWT08 . This will translate the constraint on $c_{14}$ given by
Eq.(2.14) to $0\lesssim c_{14}\lesssim\times 2.5\times 10^{-6}$, as one can
see clearly from Eq.(3.12) given in OMW18 .,

$\displaystyle 0\lesssim c_{14}\lesssim 2.5\times 10^{-5},$ (2.14)
$\displaystyle c_{14}\lesssim c_{2}\lesssim 0.095,$ (2.15)
$\displaystyle\left|c_{13}\right|\lesssim 10^{-15}.$ (2.16)

Finally, the stability of the odd-parity perturbations of BHs further requires
$c_{4}=0$ [cf. Eq.(1.1)], provided that the condition $c_{13}=0$ holds
precisely TZZW21 .

###  II.1 Particle and Universal Horizons

In this subsection, we shall provide a brief review on several horizons,
including universal horizons, as they are important to the investigations of
the QNMs of black holes in gravitational theories with broken LI. Although
they will be reviewed in the framework of the Einstein-aether theory, their
generalizations to other theories are straightforward. For more details, we
refer readers to Wang17 .

####  II.1.1 Particle Horizons

In particular, the $\ae$-theory possesses three different modes, and all of
them are moving in different speeds. In general each of these different modes
defines a particle horizon Jacobson . These horizons are the null surfaces of
the effective metrics,

$g_{\alpha\beta}^{(A)}\equiv
g_{\alpha\beta}-\left(c_{A}^{2}-1\right)u_{\alpha}u_{\beta},$ (2.17)

where $A=S,V,T$. The null surfaces for $A=S,\;V,\;T$ are called spin-0 horizon
(S0H), spin-1 horizon (S1H) and spin-2 horizon (S2H), respectively. These
three different horizons will be referred as particle horizons for the
corresponding gravitons.

In the spherically symmetric case, the static spacetimes are described by the
general metric

$\displaystyle ds^{2}$ $\displaystyle=$
$\displaystyle-F(r)dv^{2}+2B(r)dvdr+r^{2}d\Omega^{2},$ (2.18)

where $d\Omega^{2}\equiv d\theta^{2}+\sin^{2}\theta d\varphi^{2}$ and
${x}^{\mu}=(v,r,\theta,\varphi)$ denote the ingoing EF coordinates. In this
system of coordinates, the aether field takes the form

$\displaystyle u^{\alpha}\partial_{\alpha}$ $\displaystyle=$ $\displaystyle
A(r)\partial_{v}-\frac{1-F(r)A^{2}(r)}{2B(r)A(r)}\partial_{r},$ (2.19)

where $A(r)$ is an arbitrary function of $r$, and will be determined by the
Einstein-aether field equations, together with the metric coefficients $F(r)$
and $B(r)$.

Then, the out-going normal vector to the hypersurface $r=r_{0}$ is given by

$\displaystyle N_{\alpha}\equiv\frac{\partial(r-r_{0})}{\partial
x^{\alpha}}=\delta^{r}_{\alpha},$ (2.20)

where $r_{0}$ is a constant. In terms of $N_{\alpha}$, the spin-0, spin-1, and
spin-2 horizons are defined, respectively, by

$\displaystyle\left.g_{\alpha\beta}^{(A)}N^{\alpha}N^{\beta}\right|_{r=r_{A}}=0,$
(2.21)

where $r_{A}$ denotes the location of the horizon of the particle with spin-A.

On the other hand, the MH is the null surface of metric $g_{\alpha\beta}$, or
a particle horizon of $g_{\alpha\beta}^{(A)}$ with $c_{A}=1$, given by

$\displaystyle\left.g_{\alpha\beta}N^{\alpha}N^{\beta}\right|_{r=r_{\text{MH}}}=0.$
(2.22)

####  II.1.2 Universal Horizons

To define the universal horizon, let us first note that, if a BH is defined to
be a region that traps all possible causal influences, it must be bounded by a
horizon corresponding to the fastest speed. In theories with the broken LI,
the dispersion relation of a massive particle contains generically high-order
momentum terms Wang17 ,

$E^{2}=m^{2}+p_{k}^{2}k^{2}\left(1+\sum^{2(z-1)}_{n=1}{q_{n}\left(\frac{k}{M_{*}}\right)^{n}}\right),$
(2.23)

from which we can see that both of the group and phase velocities become
unbounded as $k\rightarrow\infty$, where $E$ and $k$ are the energy and
momentum of the particle considered, and $p_{k}$ and $q_{n}$’s are
coefficients, depending on the species of the particle, while $M_{*}$ is the
suppression energy scale of the higher-dimensional operators. Therefore, in
theories with the broken LI, a BH should be defined to be a region that traps
all possible causal influences, including particles with arbitrarily large
speeds $(c_{A}\rightarrow\infty)$.

Does such a region exist? To answer this question, let us first note that the
causal structure of spacetimes in such theories is quite different from that
given in GR, where the light cone at a given point $p$ plays a fundamental
role in determining the causal relationship of $p$ to other events GLLSW . In
a ultraviolet (UV) complete theory, the above dispersion relationship is
expected even in the gravitational sector 444One of such examples is the
healthy extension BPSa ; BPSb  of Hořava gravity Horava ; Wang17 , a possible
UV extension of the khronometric theory.. In such theories, the causal
structure is dramatically changed. For example, in Newtonian theory, time is
absolute and the speeds of signals are not limited. Then, the causal structure
of a given point $p$ is uniquely determined by the time difference,
$\Delta{t}\equiv t_{p}-t_{q}$, between the two events [cf. Fig. 1]. In
particular, if $\Delta{t}>0$, the event $q$ is to the past of $p$; if
$\Delta{t}<0$, it is to the future; and if $\Delta{t}=0$, the two events are
simultaneous.

Figure 1: Illustration of causal structures of spacetimes in different
theories of gravity GLLSW : (a) The light cone of the event $p$ in special
relativity. (b) The causal structure of the event $p$ in Newtonian theory.
Figure 2: Illustration of the bending of the $\phi$ = constant surfaces, and
the existence of the UH in a spherically symmetric static spacetime, where
$\phi$ denotes the globally timelike scalar field, and $t$ is the Painlevé-
Gullstrand-like coordinates, which covers the whole spacetime LSW16 .
Particles move always along the increasing direction of $\phi$. The Killing
vector $\zeta^{\mu}=\delta^{\mu}_{v}$ always points upward at each point of
the plane. The vertical dashed line is the location of the metric (Killing)
horizon, $r=r_{\text{MH}}$. The UH, denoted by the vertical solid line, is
located at $r=r_{\text{UH}}$, which is always inside the MH.

In theories with broken LI, a similar situation occurs. Thus, to consider the
causal structure of spacetimes in such theories, a globally time-like
“coordinate” needs to be introduced BS11  (See also BJS11 ). In particular,
for a given spacetime, let us first introduce a globally timelike scalar field
$\hat{\phi}$ LACW , so that

$g^{\alpha\beta}\hat{\phi}_{,\alpha}\hat{\phi}_{,\beta}>0,$ (2.24)

over the whole spacetime. It is interesting to note that, in spherically
symmetric spacetimes, the timelike aether field $u_{\mu}$ is hypersurface-
orthogonal Jacobson10 ; Wang17 , and can always be expressed as the gradient
of a scalar field Wald94 ,

$u_{\mu}=\frac{\phi_{,\mu}}{\sqrt{-\phi_{,\alpha}\phi^{,\alpha}}}.$ (2.25)

Therefore, now a natural choice is to identify $\hat{\phi}$ with $\phi$,
$\hat{\phi}=\phi$. Then, similar to Newtonian theory, this field defines
globally an absolute time direction: all particles are assumed to move along
the increasing direction of the timelike scalar field $\phi$. It is clear with
this definition of the time direction, the causality now is well-defined,
quite similar to Newtonian theory [cf. Fig. 2. See also Enrico2013 ]. With
such a definition, one can see that, for a given spacetime there may exist a
surface $r=r_{\text{UH}}$, at which the aether field $u_{\mu}$ is orthogonal
to the timelike Killing vector LGSW15  $\zeta\;(\equiv\partial_{v})$,

$\left.\zeta\cdot
u\right|_{r=r_{\text{UH}}}=-\left.\frac{1}{2A}\left(1+J\right)\right|_{r=r_{\text{UH}}}=0,$
(2.26)

where $J\equiv FA^{2}$ Chao2020b . Given that all particles move along the
increasing direction of the aether field, it is clear that a particle must
cross this surface and move inward, once it arrives at it, no matter how large
its speed is 555 Particles even with infinitely large speeds would just move
on these boundaries and cannot escape to infinity.. This is a one-way
membrane, and particles even with infinitely large speeds cannot escape from
it, once they are inside it. So, it acts as an absolute horizon to all
particles (with any speed), which is often called the UH BS11 ; BJS11 ; Wang17
.

The extension of the definition of UHs to dynamical spacetimes is given in
BMSW18 , in which it was shown explicitly that such dynamical UHs can be
formed from gravitational collapse of realistic matter fields.

###  II.2 The de Sitter-Schwarzschild Black Hole Solutions

When $c_{13}=c_{14}=0$, there exists a particular class of analytic solutions
of the background Einstein-aether field equations Chao2020b ,

$\displaystyle F(r)$ $\displaystyle=$ $\displaystyle
F_{2}\left(1-\frac{2m}{r}\right)+\frac{\Lambda}{3}r^{2},$ $\displaystyle B$
$\displaystyle=$
$\displaystyle\sqrt{F_{2}},\quad\Lambda\equiv\frac{9}{8}c_{2}w_{1}^{2}w_{2}^{2},$
$\displaystyle A(r)$ $\displaystyle=$
$\displaystyle-\frac{w_{2}}{2F}\Bigg{[}\left(\frac{1}{r^{2}}+w_{1}r\right)$
(2.27)
$\displaystyle\left.\pm\sqrt{\frac{4F}{w_{2}^{2}}+\left(\frac{1}{r^{2}}+w_{1}r\right)^{2}}\right],$

where $F_{2},\;m,\;w_{1}$ and $w_{2}$ are four integration constants. Using
the gauge freedom of the EF coordinates, without loss of the generality, we
can always set $F_{2}=1$. Then, the spacetime is precisely the anti-de Sitter
or de Sitter Schwarzschild black hole solutions, depending on the signs of the
coupling constant $c_{2}$. However, the constraints (2.14) and (2.15) tell us
that $c_{2}\geq 0$, so it must be the de Sitter Schwarzschild solution, where
$r_{s}=2m$ is the Schwarzschild radius. In this paper, we shall consider only
the case $w_{1}=0$, so that the spacetime is described precisely by the
Schwarzschild solution. Then, as $r\rightarrow\infty$, we find that

$\lim_{r\to\infty}A(r)=\mp 1,\;(F_{2}=1,\;w_{1}=0).$ (2.28)

Thus, we shall choose the “-” sign for the expression of $A(r)$ in Eq.(II.2),
so that $\lim_{r\to\infty}A(r)=+1$, which is also consistent with the
asymptotical-flatness conditions adopted in Chao2020b . Then, for a UH to
exist, we must set $w_{2}=3\sqrt{3}r_{s}^{2}/8$ LGSW15 , for which the UH is
located at

$\displaystyle r_{\text{UH}}$ $\displaystyle=$
$\displaystyle\frac{3}{4}r_{s}.$ (2.29)

Note that in the current case we have

$c_{S}=\infty,\quad c_{T}=c_{V}=1,\;(c_{13}=c_{14}=0),$ (2.30)

as can be seen from Eq. (II), which implies that now the S0H coincides with
the UH, while the S1H and S2H coincide with the MH.

It is interesting to note that in this particular case, the aether field has
no influence on the spacetime geometry. However, the aether field is non-
trivial, and quite different from the one given in the Minkowski spacetime.
This will in turn affect the perturbation equations of the Einstein-aether
gravity, as to be seen clearly in the next section.

##  III Odd-Parity Perturbations and Field Equations

To study the linear perturbations of the spherically symmetric spacetimes in
the framework of $\ae$-theory, let us first introduce the more familiar
Schwarzschild coordinate $t$ via the relation,

$\displaystyle t=v-f(r),\quad df\equiv\frac{B}{F}dr,$ (3.1)

from which we find that the background metric (2.18) and background aether
field (2.19) take the form TZZW21 ,

$\displaystyle d\bar{s}^{2}$ $\displaystyle=$
$\displaystyle-F(r)dt^{2}+\frac{B^{2}(r)}{F(r)}dr^{2}+r^{2}d\Omega^{2},$ (3.2)
$\displaystyle\bar{u}^{\alpha}\partial_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{FA^{2}+1}{2AF}\partial_{t}+\frac{FA^{2}-1}{2AB}\partial_{r}.$
(3.3)

Clearly, this metric becomes singular at the MH, $F(r_{s})=0$, so that the
coordinates cannot cover the whole region $r\in(0,\infty)$. However, now the
inner boundary is not at the MH, but rather at the UH, which is always inside
the MH. Therefore, to impose the boundary conditions at $r=r_{\text{UH}}$ and
$r=\infty$ simultaneously, we shall work out the linear perturbations in the
EF coordinates, in which the metric takes the form of Eq.(2.18). On the other
hand, setting

$\displaystyle u\equiv t-f(r),$ (3.4)

we find that the metric (3.2) takes the form

$\displaystyle d{\bar{s}}^{2}$ $\displaystyle=$
$\displaystyle-F(r)du^{2}-2B(r)dudr+r^{2}d\Omega^{2},$ (3.5)

which is the metric written in terms of the out-going EF coordinate $u$.
Although the roles of $t$ and $r$ will be exchanged across a MH, the “in-
going” and “out-going” coordinates $v$ and $u$ remain the same in the whole
spacetime. So, in the rest of the paper we always refer the “in-going” and
“out-going” in terms of the coordinates $v$ and $u$.

Once the above issues are clarified, now let us consider the perturbations of
the background solutions given by Eq.(II.2) with

$F_{2}=1,\quad w_{1}=0,\quad w_{2}=\frac{3\sqrt{3}r_{s}^{2}}{8}.$ (3.6)

For such choices, the background spacetime is precisely the Schwarzschild
vacuum solution,

$F=1-\frac{r_{s}}{r},\quad B=1,$ (3.7)

as shown in the last section.

Denoting the background metric and aether field by $\bar{g}_{\mu\nu}$ and
$\bar{u}_{\mu}$, respectively, the total metric and aether field are given by,

$\displaystyle g_{\mu\nu}$ $\displaystyle=$
$\displaystyle\bar{g}_{\mu\nu}+\epsilon h_{\mu\nu},\quad
u_{\mu}=\bar{u}_{\mu}+\epsilon w_{\mu},$ (3.8)

where $\epsilon$ is a book-marker, and we shall expand the perturbations only
to its first order. Later, we can safely set it to one. Then, working in the
in-going EF coordinates $(v,r,\theta,\phi)$, and with only the the odd-parity
part, we find that the linear perturbations $h_{\mu\nu}$ and $w_{\mu}$ can be
cast in the forms Thomp2017 ; TZZW21 ,

$\displaystyle h_{\mu\nu}$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\begin{pmatrix}0&0&C_{lm}\csc\theta\partial_{\varphi}&-C_{lm}\sin\theta\partial_{\theta}\\\
0&0&J_{lm}\csc\theta\partial_{\varphi}&-J_{lm}\sin\theta\partial_{\theta}\\\
sym&sym&G_{lm}\csc\theta\big{(}\cot\theta\partial_{\varphi}-\partial_{\theta}\partial_{\varphi}\big{)}&sym\\\
sym&sym&\frac{1}{2}G_{lm}\big{(}\sin\theta\partial^{2}_{\theta}-\cos\theta\partial_{\theta}-\csc\theta\partial^{2}_{\varphi}\big{)}&-G_{lm}\sin\theta\big{(}\cot\theta\partial_{\phi}-\partial_{\theta}\partial_{\varphi}\big{)}\end{pmatrix}Y_{lm}(\theta,\varphi),$

and

$\displaystyle w_{\mu}$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\begin{pmatrix}0\\\ 0\\\
a_{lm}\csc\theta\partial_{\varphi}\\\
-a_{lm}\sin\theta\partial_{\theta}\end{pmatrix}Y_{lm}(\theta,\varphi),$ (3.10)

where $Y_{lm}(\theta,\varphi)$ stands for the spherical harmonics, and
$C_{lm},\;G_{lm},\;J_{lm}$ and $a_{lm}$ are functions of $v$ and $r$ only 666
It must not be confused with the function $J_{lm}$ introduced here and the
tensor $J_{\alpha\beta}$ appearing in Eq.(II).. Note that when calculating the
field equations we will set $m=0$ in the above expressions so that
$\partial_{\varphi}Y_{lm}(\theta,\varphi)=0$, as now the background has the
spherical symmetry, and the corresponding linear perturbations do not depend
on $m$ Regge57 ; Thomp2017  777Notice that, since we are using the EF
coordinate and following some different conventions, the terms
$\\{C_{lm},J_{lm},G_{lm}\\}$ in (III) are not necessarily equal to their
counterparts given in Thomp2017 ..

###  III.1 Gauge Transformations

For later convenience, we first investigate the infinitesimal gauge
transformations (Recall that we only consider the odd-parity perturbations and
$m=0$)

$\displaystyle x^{\alpha}\to
x^{\prime\alpha}=x^{\alpha}+\epsilon\xi^{\alpha},$ (3.11)

where

$\displaystyle\xi^{\alpha}$ $\displaystyle=$
$\displaystyle-\frac{\csc\theta\partial_{\theta}Y_{lm}(\theta,\phi)}{r^{2}}\left\\{0,0,0,1\right\\}\xi,~{}~{}~{}~{}$
(3.12)

with $\xi$ being a function of $v$ and $r$. Under the transformation of
Eq.(3.11), we have

$\displaystyle\Delta
w_{\mu}\equiv\left(w_{\mu}\right)_{new}-\left(w_{\mu}\right)_{old}=-{\cal{L}}_{\xi}{\bar{u}}_{\mu},$
$\displaystyle\Delta
h_{\mu\nu}\equiv\left(h_{\mu\nu}\right)_{new}-\left(h_{\mu\nu}\right)_{old}=-{\cal{L}}_{\xi}{\bar{g}}_{\mu\nu},$
(3.13)

where $\cal{L}$ stands for the Lie derivative Invb . From Eq.(III.1) we find
TZZW21

$\displaystyle\Delta
C_{lm}\equiv\left(C_{lm}\right)_{old}-\left(C_{lm}\right)_{new}={\dot{\xi}},$
$\displaystyle\Delta
G_{lm}\equiv\left(G_{lm}\right)_{old}-\left(G_{lm}\right)_{new}=-2\xi,$
$\displaystyle\Delta
J_{lm}\equiv\left(J_{lm}\right)_{old}-\left(J_{lm}\right)_{new}=\xi^{\prime}-\frac{2}{r}\xi,$
$\displaystyle\Delta
a_{lm}\equiv\left(a_{lm}\right)_{old}-\left(a_{lm}\right)_{new}=0,$ (3.14)

where a prime and a dot stand for the derivatives with respect to $r$ and $v$,
respectively. With the above gauge transformations, we can construct the
gauge-invariant quantities, and due to the presence of the aether field, three
such independent quantities can be constructed, in contrast to the
relativistic case, in which only two such quantities can be constructed. These
three gauge invariants can be defined as

$\displaystyle{\cal{X}}_{lm}(v,r)$ $\displaystyle\equiv$
$\displaystyle{C}_{lm}+\frac{1}{2}{\dot{G}}_{lm},$
$\displaystyle{\cal{Y}}_{lm}(v,r)$ $\displaystyle\equiv$
$\displaystyle\frac{2}{r}C_{lm}-C_{lm}^{\prime}+{\dot{J}}_{lm},$
$\displaystyle{\cal{Z}}_{lm}(v,r)$ $\displaystyle\equiv$ $\displaystyle
a_{lm}.$ (3.15)

Of course, any combination of these quantities is also gauge-invariant. By
properly fixing the gauge, the resultant field equations will be simplified
considerably. In the current case, we find that one of the most convenient
gauge choices is

$C_{lm}=0,$ (3.16)

from which we find

$\displaystyle\xi(v,r)$ $\displaystyle=$
$\displaystyle\int\left[C_{lm}(v,r)\right]_{old}dv+{\hat{C}}_{1}(r),$ (3.17)

where ${\hat{C}}_{1}$ is an arbitrary function of $r$. Therefore, with this
choice of gauge, we still have the gauge residual

${\hat{\xi}}(v,r)={\hat{C}}_{1}(r).$ (3.18)

We shall come back to this point later, and fix the gauge residual completely
by properly imposing an additional gauge condition.

###  III.2 Linearized Field Equations

When the spacetimes are in the vacuum, we have $T_{\mu\nu}=0,\;T_{\mu}=0$, and
then the field equations (2.7) and (2.8) reduce to

$\displaystyle E_{\mu\nu}\equiv G_{\mu\nu}-S_{\mu\nu}=0,$ (3.19)
$\displaystyle\AE^{\mu}=0,$ (3.20)

where $G_{\mu\nu}\equiv R_{\mu\nu}-g_{\mu\nu}R/2$. To the first-order of
$\epsilon$ and with the gauge (3.16), we find that there are only four non-
trivial equations, which are given by

$E_{\phi t}=E_{\phi r}=E_{\phi\theta}=\AE^{\phi}=0,$ (3.21)

and can be cast respectively in the forms given by Eqs.(A. 2) - (A. 5) in
Appendix A.

It is interesting to note that the coefficients of these differential
equations, i.e., $\rho_{abc}$ given by Eqs.(Appendix A: The coefficients of
$\rho_{abc}$, $\kappa_{n}$ and $\alpha_{m}$) - (Appendix A: The coefficients
of $\rho_{abc}$, $\kappa_{n}$ and $\alpha_{m}$), are functions of $r$, $c_{1}$
and $l$ only, and some of them contain no $c_{1}$. More importantly, from Eq.
(A. 5) we can see clearly that $a_{lm}$ has been decoupled from $J_{lm}$ and
$G_{lm}$.

###  III.3 Master Equations

Since there are only three independent unknowns in Eqs. (A. 2)-(A. 5),
$G_{lm},\;J_{lm}$ and $a_{lm}$, one of the above four equations must not be
independent, although this non-dependent equation can still be used to
simplify the derivations of the master equation. To simplify the following
discussions, using the rescaling symmetry of the background Chao2020b , we
first set $r_{s}=1$. Then, we have

$f(r)\equiv\int{\frac{B}{F}dr}=r+\ln\left|r-1\right|.$ (3.22)

To find the master equations for the gauge invariants introduced in
Eq.(III.1), we first eliminate the ${\ddot{a}}_{lm}$ terms from Eqs.(A. 2) and
(Appendix A: The coefficients of $\rho_{abc}$, $\kappa_{n}$ and $\alpha_{m}$).
Then, combining this resultant equation with Eq.(A. 4) we find that the terms
proportional to ${\dot{G}}_{lm}^{\prime}$ and ${\dot{G}}_{lm}$ can be made
disappear. Combining this resultant equation with Eq.(A. 2) and the first-
order $v$-derivative of Eq.(Appendix A: The coefficients of $\rho_{abc}$,
$\kappa_{n}$ and $\alpha_{m}$), we obtain an equation that contains the
${\dddot{a}}_{lm}$ and ${a}_{lm}^{\prime\prime}$ terms. Finally, combining
this new resultant equation with Eq.(A. 5) and its $v$ derivative, we arrive
at the master equation

$\displaystyle\Bigg{\\{}\frac{(r-1)^{2}}{r^{2}}\frac{\partial^{2}}{\partial
r^{2}}+\frac{(r-1)(2r-1)}{r^{3}}\frac{\partial}{\partial
r}+2\left(1-\frac{1}{r}\right)\frac{\partial^{2}}{\partial r\partial
v}+\frac{2(r-1)}{r^{2}}\frac{\partial}{\partial
v}-\frac{(r-1)[l(l+1)r-4]}{r^{4}}\Bigg{\\}}{\cal{Y}}_{lm}=0,~{}~{}~{}~{}$
(3.23)

where ${\cal{Y}}_{lm}=\dot{J}_{lm}$. Setting

$\displaystyle{\cal{Y}}_{lm}(v,r)$ $\displaystyle=$ $\displaystyle e^{-i\omega
v}{\bar{\cal{Y}}}_{lm}(r),$ (3.24)

where $\omega$ is the QNM frequency to be calculated below, (3.23) reduces to

$\displaystyle\Bigg{\\{}\frac{(r-1)^{2}}{r^{2}}\frac{\partial^{2}}{\partial
r^{2}}+\left[\frac{(r-1)(2r-1)}{r^{3}}-2i\omega\left(1-\frac{1}{r}\right)\right]\frac{\partial}{\partial
r}-\frac{2i\omega(r-1)}{r^{2}}-\frac{(r-1)[l(l+1)r-4]}{r^{4}}\Bigg{\\}}\bar{\cal{Y}}_{lm}=0.~{}~{}~{}~{}$
(3.25)

To eliminate the first-order derivative term of the above equation, following
the general prescription presented in Appendix B, we introduce the function
$\Psi(r)$ via the relation,

$\displaystyle\Psi(r)$ $\displaystyle\equiv$ $\displaystyle e^{-i\omega
f(r)}r{\bar{\cal{Y}}}_{lm}(r),$ (3.26)

so that Eq. (3.25) can be cast in the simple form,

$\displaystyle\left[\frac{d^{2}}{dx^{2}}+\left(\omega^{2}-V_{g}\right)\right]\Psi=0,$
(3.27)

where $x\equiv f(r)$ given by Eq.(3.22), and

$\displaystyle V_{g}$ $\displaystyle\equiv$
$\displaystyle\frac{(r-1)[l(l+1)r-3]}{r^{4}}.$ (3.28)

Note that, the range of the new variable $x\in(-\infty,+\infty)$ covers only
the region $r\in(1,\infty)$. Thus, as will be seen later, we have to look for
another variable to cover the whole region $r\in\left(r_{UH},\infty\right)$
where $r_{UH}=3/4$. However, this equation is very useful when we study the
asymptotical behavior of the perturbations at the spatial infinity
$r\rightarrow\infty$, as can be seen from the discussions to be presented
below in Section IV.A. It is interesting to note that $V_{g}$ is identical to
the effective potential given in the Regge-Wheeler equation in GR Kono2011 .
However, due to the different locations of the inner boundaries, in general
different BH spectroscopies are expected.

On the other hand, solving for ${\dot{G}}_{lm}$ and $a_{lm}^{\prime\prime}$
from Eqs.(A. 2) and (A. 5) we find

$\displaystyle{\cal{X}}_{lm}$ $\displaystyle=$
$\displaystyle-\frac{1}{(l+2)(l-1)}$ (3.29)
$\displaystyle\times\left[2(r-1)+r(r-1)\frac{\partial}{\partial
r}+r^{2}\frac{\partial}{\partial
v}\right]{\cal{Y}}_{lm},~{}~{}~{}~{}~{}~{}~{}$

where Eq.(III.1) was used in writing down the above expression. Then, from
Eqs.(III.1) and (3.16) we find

$\displaystyle G_{lm}(v,r)$ $\displaystyle=$ $\displaystyle
2\int{{\cal{X}}(v,r)dv}+\hat{G}_{lm}(r),$ (3.30)

where $\hat{G}_{lm}(r)$ is an arbitrary function of $r$ only. However, using
the gauge residual (3.18), we can always set it to zero, by choosing

$\hat{C}(r)=-\frac{1}{2}\hat{G}_{lm}(r).$ (3.31)

Clearly, with this choice the remaining gauge freedom is completely fixed.

In addition, following Appendix B we find that Eq.(A. 5) can be written in the
form

$\displaystyle\left[\frac{(r-1)^{2}}{r^{2}}\frac{d^{2}}{dr^{2}}+\frac{r-1}{r^{3}}\frac{d}{dr}-V_{\text{eff}}(r)\right]\psi$
$\displaystyle=$ $\displaystyle 0,~{}~{}~{}~{}$ (3.32)

where

$\displaystyle\psi(r)$ $\displaystyle\equiv$
$\displaystyle\sqrt{\frac{r}{r-1}}\exp\left\\{\int\frac{54(1-r)+r^{2}\left[128(r-1)r+3i\sqrt{768(r-1)r^{3}+81}\omega_{\ae}\right]}{(3-4r)^{2}(r-1)r[8r(2r+1)+3]}dr\right\\}{\hat{\cal{Z}}}_{lm}(r),$
(3.33) $\displaystyle V_{\text{eff}}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{4\left(256r^{8}-256r^{7}+27r^{4}\right)}\Big{[}1024l(l+1)r^{6}-1024\left(2l^{2}+2l+1\right)r^{5}+256\left(4l^{2}+4l+9\right)r^{4}$
(3.34)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-1280r^{3}+2160r^{2}-4212r+2079\Big{]},$

with

$\displaystyle a_{lm}(v,r)=e^{-i\omega_{\ae}t}{\hat{\cal{Z}}}_{lm}(r).$ (3.35)

Here $t$ is related to $v$ through Eqs.(3.1) and (3.22). Notice that, a
subscript $\ae$ is added to $\omega$ here to allow it to be different from the
QNM spectra $\omega$ obtained from Eq.(3.27).

It is remarkable to note that $a_{lm}$ depends on $\omega_{\ae}$ only through
Eq.(3.33), while Eq.(3.32) shows that $\psi$ itself does not depend on
$\omega_{\ae}$, although its boundary conditions will depend on $\omega_{\ae}$
through $a_{lm}$. We shall come back to this point in the next section.

##  IV QNMs of black holes in Einstein-aether Theory

From the last section, we can see that solving the linearized Einstein-aether
field equations for the odd-parity perturbations now reduces to solving the
master equations (3.27) and (3.32) with proper boundary conditions. As in GR,
they have solutions only for some particular choices of $\omega$ and
$\omega_{\ae}$, which form the spectra of the QNMs of the corresponding BH.

In contrast to GR, due to the existence of superluminal propagations caused by
the LI violation, the location of the inner boundary conditions now is moved
from MHs to UHs, which are always inside the MHs. There are at least two
aspects that make us to believe that QNMs in $\ae$-theory in general deviate
from those of GR:

  * •

The different locations of the inner boundaries. It is true that in both
theories it is all required that only “in-going” waves be allowed at the inner
boundaries. However, when these conditions are imposed at UHs (as required by
$\ae$-theory), it is expected that at the corresponding MHs both “in-going”
and “out-going” waves exist (which is quite different from that of GR),
because of the existence of the superluminal propagations.

  * •

The different dynamical differential equations. Despite the fact that the
background in both theories is the same and all described by the Schwarzschild
BH, the linearized perturbation equations are different, and in general it is
expected that these also alternate the BH QNMs.

With the above in mind, let us turn to consider the initial conditions
respectively at the UH ($r=r_{\text{UH}}=3/4$) and spatial infinity
($r\to+\infty$). Once these are done, we turn ourselves to consider the
smoothness conditions at the MH.

###  IV.1 Boundary Conditions at the Spatial Infinities

As mentioned above, the boundary conditions at the spatial infinity are that
only “out-going” waves exist. To see the implication of these conditions on
the gauge-invariant quantities ${\cal{Y}}_{lm}$ and $a_{lm}$, let us first
note that as $r\to+\infty$, to the leading order, we have $f(r)\simeq
r,\;V_{g}(r)\rightarrow 0$, and $v\simeq t+r$ and $u\simeq t-r$. Then,
Eq.(3.27) has the general solutions,

$\Psi|_{r\to\infty}=a_{+}e^{i\omega x}+a_{-}e^{-i\omega x},$ (4.1)

where $a_{\pm}$ are two integration constants. Inserting the two branches of
(4.1) into (3.26), and together with Eq. (3.24), we find that

$\displaystyle\left.{\cal{Y}}_{lm}(v,r)\right|_{r\to+\infty}$
$\displaystyle\propto$ $\displaystyle ra_{\pm}e^{-i\omega(t\mp x)}.$ (4.2)

Therefore, to have only outgoing-waves, we must set $a_{-}=0$. As a result,
when $r$ is very large, we expect $\Psi(r)$ to take the form,

$\displaystyle\Psi$ $\displaystyle=$ $\displaystyle e^{i\omega
x}\sum_{n=0}^{\infty}\frac{a_{n}}{r^{n}},$ (4.3)

where $a_{n}$’s are constants to be determined with the Frobenius-like method
Butkov . In particular, inserting the above expression into Eq.(3.27), we find
the following recursion relation,

$\displaystyle 2in\omega a_{n}-\Big{[}(n-1)(n+2i\omega)-l^{2}-i\Big{]}a_{n-1}$
$\displaystyle-\Big{(}l^{2}+l-2n^{2}+5n+1\Big{)}a_{n-2}+(n-4)na_{n-3}=0,$

from which we can write all $a_{n}$’s ($n\geq 1$) in terms of $a_{0}$ with
$a_{n}\propto a_{0}$. As will be seen later, the amplitude of $\Psi(r)$ has no
contributions to our main results, Therefore, without loss of the generality,
we can always set $a_{0}=1$.

###  IV.2 Boundary Conditions at Universal Horizons

On the other hand, as $r\to r_{\text{UH}}\;(=3/4)$, we note that Eq. (3.25)
can be written as

$\displaystyle\left[{\tilde{W}}_{g}(r)\frac{d^{2}}{dy^{2}}+\left(\omega^{2}-{\tilde{V}}_{g}\right)\right]{\tilde{\Psi}}=0,$
(4.5)
$\displaystyle{\tilde{\Psi}}\equiv\sqrt{\frac{r(1-r)}{{\tilde{p}}}}e^{-i\omega
f}{\bar{\cal{Y}}}_{lm},\quad{\tilde{W}}_{g}\equiv\frac{\eta_{1}}{{\tilde{p}}^{2}},$
(4.6)

where

$\displaystyle\frac{dr}{dy}$ $\displaystyle\equiv$
$\displaystyle{\tilde{p}}(r)\equiv\frac{\sqrt{l^{2}+l-4}+12+l(l+1)(4r-3)-16r}{3\sqrt{l^{2}+l-4}+12+l(l+1)(4r-3)-16r},$
$\displaystyle{\tilde{V}}_{g}$ $\displaystyle\equiv$
$\displaystyle-\left\\{\frac{{\tilde{W}}_{g}}{4}\left[2{\tilde{p}}\frac{d\left({\tilde{p}}^{\prime}-{\tilde{p}}\eta_{2}/\eta_{1}\right)}{dr}-\left({\tilde{p}}^{\prime}-\frac{{\tilde{p}}\eta_{2}}{\eta_{1}}\right)^{2}\right]+\eta_{3}\right\\},$
$\displaystyle\eta_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{(r-1)^{2}}{r^{2}},\quad\eta_{2}\equiv\frac{(r-1)(2r-1)}{r^{3}},\quad\eta_{3}\equiv-\frac{(r-1)\left(l^{2}r+lr-4\right)}{r^{4}}.$
(4.7)

Thus, as $r\to r_{\text{UH}}$, we have

$\displaystyle\frac{d^{2}{\tilde{\Psi}}}{dy^{2}}+\omega^{2}{\tilde{\Psi}}$
$\displaystyle=$ $\displaystyle 0,$ (4.8)

which leads to

$\displaystyle{\tilde{\Psi}}|_{r\to r_{\text{UH}}}$ $\displaystyle=$
$\displaystyle b_{+}e^{+i\omega y}+b_{-}e^{-i\omega y},$ (4.9)

with $b_{\pm}$ being arbitrary constants. Combining Eqs. (4.9), (4.6) with Eq.
(3.24), we find that

$\displaystyle\left.{\cal{Y}}_{lm}(t,r)\right|_{r\to r_{\text{UH}}}$
$\displaystyle\propto$ $\displaystyle b_{\pm}e^{-i\omega\left(t\mp
y\right)}\sqrt{\frac{{\tilde{p}}}{r(1-r)}}$ (4.10) $\displaystyle\propto$
$\displaystyle b_{\pm}e^{-i\omega\left(t\mp r\right)}.$

Therefore, to satisfy the pure ingoing-wave conditions, we must set $b_{+}=0$.
As a result, when $r\to r_{\text{UH}}$, we expect ${\tilde{\Psi}}(r)$ to take
the form,

$\displaystyle{\tilde{\Psi}}$ $\displaystyle=$ $\displaystyle e^{-i\omega
y}\sum_{n=0}^{\infty}b_{n}\left(r-\frac{3}{4}\right)^{n}.$ (4.11)

Similar to $a_{n}$’s, inserting the above expression into Eq.(4.5) we shall
find a recursion relation for $b_{n}$’s, and from which we can write all
$b_{n\geq 1}$’s in terms of $b_{0}$ with $b_{n}\propto b_{0}$. Again, as will
be seen, the choice of $b_{0}$ is actually irrelevant to our solutions of
QNMs, so without loss of the generality, we can also set it to one, $b_{0}=1$.

###  IV.3 Smoothness Conditions across Metric Horizons

In principle, once the two boundary conditions are given, we can solve the two
master equations (3.27) and (3.32) to find out the spectrum of $\omega$
($\omega_{\ae}$) for any given $l$. However, since both of them are singular
at the MH, we have to solve them first in the regions
$r\in(r_{UH},1-\epsilon)$ and $r\in(1+\epsilon,\infty)$ separately, and then
match the solutions together across the MH, where $0<\epsilon\ll 1$ [Notice
that $\epsilon$ used here must not confused with the ones appearing in
Eq.(3.8)].

To find out the proper matching conditions in the neighborhood of $r=1$,
Eq.(3.27) cannot be used, as it is valid only for $r\in(1,\infty)$. To
overcome this problem, let us first replace ${\bar{\cal{Y}}}_{lm}$ appearing
in Eq. (3.25) by

${\hat{\cal{Y}}}_{lm}\equiv e^{-i\omega f}{\bar{\cal{Y}}}_{lm},$ (4.12)

so that it reduces to

$\displaystyle\Bigg{\\{}\frac{(r-1)^{2}}{r^{2}}\frac{d^{2}}{dr^{2}}+\frac{(r-1)(2r-1)}{r^{3}}\frac{d}{dr}+\left[\omega^{2}-\frac{(r-1)\left(l^{2}r+lr-4\right)}{r^{4}}\right]\Bigg{\\}}{\hat{\cal{Y}}}_{lm}(r)=0.~{}~{}~{}~{}$
(4.13)

Then, we expand $\hat{\cal{Y}}_{lm}(r)$ as

$\displaystyle\hat{\cal{Y}}_{lm}(r)=(r-1)^{s}\sum_{n=0}^{\infty}d_{n}\left(r-1\right)^{n},$
(4.14)

where $d_{n}$’s are constants to be determined with the Frobenius method. In
particular, inserting the above expression into Eq.(4.13), we find two
branches of solutions, given, respectively, by

$s_{\pm}=\pm i\omega.$ (4.15)

Inserting it together with Eq.(4.14) into Eq.(4.13) we find

$\displaystyle
n(n+2i\omega)d^{+}_{n}+\left[-l^{2}-l+2n^{2}+n(-3+4i\omega)+2\omega^{2}-3i\omega+5\right]d^{+}_{n-1}$
$\displaystyle+\left[-l^{2}-l+n^{2}+n(-3+2i\omega)+5\omega^{2}-3i\omega+2\right]d^{+}_{n-2}+4\omega^{2}d^{+}_{n-3}+\omega^{2}d^{+}_{n-4}=0,$
(4.16)

for $s=s_{+}$, and

$\displaystyle
n(n-2i\omega)d^{+}_{n}+\left[-l^{2}-l+2n^{2}+n(-3-4i\omega)+2\omega^{2}+3i\omega+5\right]d^{+}_{n-1}$
(4.17)
$\displaystyle+\left[-l^{2}-l+n^{2}+n(-3-2i\omega)+5\omega^{2}+3i\omega+2\right]d^{+}_{n-2}+4\omega^{2}d^{+}_{n-3}+\omega^{2}d^{+}_{n-4}=0,$

for $s=s_{-}$. Then, similar to $a_{n}$’s given by Eq.(IV.1), all these
coefficients can be written in terms of $d^{\pm}_{0}$. With the same reasons
as given for the choice of $a_{0}$ and $b_{0}$, without loss of the
generality, we can always set $d^{\pm}_{0}=1$. Therefore, for any given
$\omega$ and $l$, we find the general solution of $\hat{\cal{Y}}_{lm}(r)$ in
the neighborhood of $r=1$,

$\displaystyle\hat{\cal{Y}}_{lm}(r)$ $\displaystyle=$
$\displaystyle(r-1)^{s_{-}}\sum_{n=0}^{\infty}d_{n}^{-}\left(r-1\right)^{n}+{\cal{A}}(r-1)^{s_{+}}\sum_{n=0}^{\infty}d_{n}^{+}\left(r-1\right)^{n},$
(4.18)

with ${\cal{A}}$ being a complex number, and $d_{n}^{\pm}$’s are given by
Eqs.(IV.3) and (4.17) with $d^{\pm}_{0}=1$.

It should be noted that in practice Eqs. (4.3), (4.11) and (4.18) cannot be
expanded to infinite orders. Instead, we shall expand them to the orders by
imposing certain tolerances to the numerical errors of our numerical
computations, and make sure such orders are high enough, so that our main
results and conclusions will remain the same, i.e., the numerical solutions
adopt convergent behaviors.

###  IV.4 Algorithm for Solving Master Equations

To solve Eq.(3.25), following Chandra1975  we introduce the mode function
$\Phi(r)$ by

$\displaystyle{\hat{\cal{Y}}}_{lm}(r)$ $\displaystyle=$
$\displaystyle\text{Exp}\left(i\int^{r}\Phi(r)dr\right),$ (4.19)

or equivalently

$\displaystyle\Phi(r)$ $\displaystyle=$
$\displaystyle-i\frac{d\ln{\hat{\cal{Y}}}_{lm}}{dr}.$ (4.20)

Then, Eq.(3.25) reduces to

$\displaystyle\frac{i(r-1)^{2}}{r^{2}}\Phi^{\prime}-\frac{(r-1)^{2}}{r^{2}}\Phi^{2}+\frac{i\left(2r^{2}-3r+1\right)}{r^{3}}\Phi$
$\displaystyle-\frac{(r-1)\left(l^{2}r+lr-4\right)}{r^{4}}+\omega^{2}=0.$
(4.21)

It is remarkable to note that such introduced $\Phi(r)$ does not depend on the
amplitudes of ${\hat{\cal{Y}}}_{lm}$ around any of the points
$r=\left\\{r_{\pm},\;r_{\text{UH}},\;r_{\text{max}}\right\\}$, with
$r_{\pm}\equiv 1\pm\epsilon$ ( $0<\epsilon\ll 1$) and $r_{\text{max}}$ a
positive number to be specified. As a result, we can assign $a_{0},\;b_{0}$
and $d_{0}^{\pm}$ any values. As mentioned previously, without loss of the
generality, we shall set all of them to one.

For our convenience, we further introduce a new variable $z$ by $z\equiv
2r/(r+3/4)$ in Eq.(IV.4). In this way, $r\in[3/4,1)\cup(1,\infty)$ is mapped
to $z\in[1,8/7)\cup(8/7,2)$. Therefore, Eq.(IV.4) becomes

$\displaystyle\frac{2i\left(7z^{2}-22z+16\right)^{2}}{27z^{2}}\frac{d\Phi}{dz}-\frac{(8-7z)^{2}}{9z^{2}}\Phi^{2}$
$\displaystyle-\frac{8i\left(35z^{3}-138z^{2}+168z-64\right)}{27z^{3}}\Phi+\omega^{2}$
$\displaystyle-\frac{16(7z-8)(z-2)^{2}\left[\left(3l^{2}+3l+16\right)z-32\right]}{81z^{4}}=0.$

Thus, we obtain a complex ordinary differential equation (ODE) of the first
order, which will be much easier to solve than Eq.(3.25). To solve Eq.(IV.4)
we use the shooting method. Specifically, for any given $l$, we integrate
Eq.(IV.4) as follows:

  1. 1.

With a good guess, we first choose a (complex) value of $\omega$. Then, by
using the polynomial solution (4.11), we create the “initial value”
${\tilde{\Psi}}_{\text{UH}}(\omega)\equiv{\tilde{\Psi}}(r=r_{\text{UH}},\omega)$.
After that, we calculate
$\Phi_{\text{UH}}(\omega)\equiv\Phi(z=z_{\text{UH}},\omega)$ using Eq.(4.20)
[and Eqs.(4.12), (4.6)], with
$z_{\text{UH}}=2r_{\text{UH}}/(r_{\text{UH}}+3/4)$.

  2. 2.

With the same $\omega$, using the polynomial solution (4.3), we create the
“initial value”
$\Psi_{\text{max}}(\omega)\equiv\Psi(r=r_{\text{max}},\omega)$. After that, we
calculate $\Phi_{\text{max}}(\omega)\equiv\Phi(z=z_{\text{max}},\omega)$ using
Eq.(4.20) [and Eqs.(4.12), (3.26)], with
$z_{\text{max}}=2r_{\text{max}}/(r_{\text{max}}+3/4)$. Here $r_{\text{max}}$
is chosen to be sufficiently large.

  3. 3.

By using the initial value $\Phi_{\text{UH}}(\omega)$, integrate Eq.(IV.4)
from $z_{\text{UH}}$ to $z_{-}=r_{-}/(r_{-}+3/4)$ to obtain
$\Phi_{\text{-}}\equiv\Phi(z_{\text{-}})$. Equalizing this $\Phi_{\text{-}}$
and its counterpart given by the polynomial (4.18) [and Eq.(4.20)], we can
solve for the combination coefficient ${\cal A}$, to be denoted by ${\cal
A}_{-}$.

  4. 4.

By using the initial value $\Phi_{\text{max}}(\omega)$, integrate Eq.(IV.4)
from $z_{\text{max}}$ to $z_{+}=r_{+}/(r_{+}+3/4)$ to obtain
$\Phi_{\text{+}}\equiv\Phi(z_{\text{+}})$. Equalizing this $\Phi_{\text{+}}$
and its counterpart given by the polynomial (4.18) [and Eq.(4.20)], we can
again solve for a combination coefficient ${\cal A}$, to be denoted by ${\cal
A}_{+}$.

  5. 5.

Logically repeat steps 1 - 4 to search for the desired $\omega$, with which we
have $\Delta{\cal A}\equiv|{\cal A}_{-}-{\cal
A}_{+}|\lesssim\epsilon_{\text{max}}$, where $\epsilon_{\text{max}}$
represents our tolerance to the uncertainty. Clearly, this $\Delta{\cal A}$
could be seen as an indicator to the accuracy of our calculations. Note that
these iterations could be executed automatically by using Mathematica.

Table 1: Usages of equations in generating the suitable boundary conditions for solving Eq.(IV.4). Type of boundary condition |  |  Characterized by |  |  Polynomial |  |  Translated to $\Phi(r)$ with
---|---|---|---|---|---|---
|  |  |  | expression |  |
Pure out-going wave at $r\to+\infty$ |  |  Eq.(4.1), with $a_{-}=0$ |  | Eq.(4.3) |  | Eqs.(4.20), (4.12), (3.26)
Pure in-going wave at $r\to(3/4)^{+}$ |  |  Eq.(4.9), with $b_{+}=0$ |  | Eq.(4.11) |  | Eqs.(4.20), (4.12), (4.6)
Smoothness across $r=1$ |  | N/A |  | Eq.(4.18) |  | Eq.(4.20)

###  IV.5 Numerical Results

Using the algorithm introduced in last subsection, we can obtain eigenvalues
of $\omega$ for any given $l$. However, since the algorithm introduced in this
section is quite algebraically involved, before presenting our results, for
the sake of readers’ convenience, let us first briefly provide a summary to
the pivot equations we use in the shooting method. Basically, our main target
is to obtain eigenvalues of $\omega$ by solving the first-order non-linear
complex differential equation (IV.4). For this purpose, suitable boundary
conditions are needed. We summarize the relative necessary equations in Table
1. Eq.(IV.4) plus the ones given in Table 1 form a complete system, and the
eigenvalues of $\omega$ given in Table 2 are found by using the shooting
method with the assistance of Mathematica. In this table we also provide the
corresponding results from GR (see, e.g., Chao2022 ). Comparing the results
from $\ae$-theory and GR, we can find the differences between these two
theories. To characterize these discrepancies, let’s define Berti2009 ;
Ghosh2021

$\displaystyle\delta f$ $\displaystyle\equiv$
$\displaystyle\frac{Re(\omega^{\text{\ae}})}{Re(\omega^{\text{GR}})}-1,$
$\displaystyle\delta\tau$ $\displaystyle\equiv$
$\displaystyle\frac{Im(\omega^{\text{GR}})}{Im(\omega^{\text{\ae}})}-1,$
(4.23)

which are also given for different modes in Tabel 2. Using the current data
provided by the LVK scientific collaboration Ghosh2021 ; LVK2022 , it is hard
to impose further constraints on $\ae$-theory. For instance, combining the
results from Ghosh2021 ; LVK2022 , we find that the most stringent constraint
on $\delta f$ is right at the order ${\cal{O}}(10^{-2})$, which is similar to
the largest discrepancy that we have presented in Table 2. Nonetheless, such
discrepancies should be well within the detectability of LISA, TianQin and
Taiji. In particular, according to Table II provided in Liu2020 , the
sensitivity of LISA and TianQin to, e.g., $\delta f$, could reach about
${\cal{O}}(10^{-5})$. Therefore, they are very likely to be able to
distinguish the results from the $\ae$-theory and GR, as listed in Table 2.

Table 2: Values of $\omega$ for any given $l$ in $\ae$-theory and GR, where we have set $c=G_{N}=r_{MH}(=r_{s})=1$. $l$ |  | $\ae$ |  | $\Delta{\cal A}$ |  | GR |  | $\delta f$ |  | $\delta\tau$
---|---|---|---|---|---|---|---|---|---|---
$l=2$ |  |  $-0.74314-0.17787i$ |  | $4.76\times 10^{-34}$ |  |  $-0.74734-0.17792i$ |  | -0.00562 |  | -0.00030
|  |  $-0.67267-0.54966i$ |  | $4.54\times 10^{-30}$ |  | $-0.69342-0.54783i$ |  | -0.02993 |  | 0.00334
|  |  $-0.56544-0.91304i$ |  | $9.38\times 10^{-27}$ |  | $-0.60211-0.95656i$ |  | -0.06090 |  | -0.04549
$l=3$ |  |  $-1.19890-0.18506i$ |  | $1.40\times 10^{-33}$ |  | $-1.19889-0.18541i$ |  | 0.00001 |  | -0.00189
|  |  $-1.16448-0.56037i$ |  | $3.69\times 10^{-30}$ |  | $-1.16529-0.56260i$ |  | -0.00070 |  | -0.00396
|  |  $-1.07995-0.95446i$ |  | $3.28\times 10^{-27}$ |  | $-1.10337-0.95819i$ |  | -0.02122 |  | -0.00388
$l=4$ |  |  $-1.61838-0.18832i$ |  | $1.05\times 10^{-18}$ |  | $-1.61836-0.18833i$ |  | 0.00002 |  | -0.00006
|  |  $-1.59342-0.56847i$ |  | $1.49\times 10^{-29}$ |  | $-1.59326-0.56867i$ |  | 0.00010 |  | -0.00034
|  |  $-1.54705-0.95738i$ |  | $9.75\times 10^{-25}$ |  | $-1.54542-0.95982i$ |  | 0.00106 |  | -0.00254
$l=10$ |  |  $-3.98329-0.57582i$ |  | $2.14\times 10^{-14}$ |  | $-3.98329-0.57582i$ |  | 0.00000 |  | 0.00000

##  V  Instability of the aether field

In the last section, we consider the metric perturbations only, represented by
the master equation (3.25). In this section, we turn to study the aether
perturbations given by Eqs.(3.32) - (3.35). To this goal, inserting Eq.(3.35)
into Eq.(A. 5) we obtain

$\displaystyle\left\\{\kappa_{1}\frac{d^{2}}{dr^{2}}+\kappa_{2}\frac{d}{dr}+\kappa_{3}\right\\}{\hat{\cal{Z}}}_{lm}$
$\displaystyle=$ $\displaystyle 0,$ (5.1)

where $\kappa_{1,2,3}$ are given in Appendix A. Note that (5.1) is at the same
position as (4.13). To solve it, following what we did in the last section, we
first specify its boundary conditions at the universal horizon and spatial
infinity, as well as the smoothness conditions across the metric horizon, and
then solve it with such conditions by the shooting method, one side from the
UH to the MH, and the other side from $r_{\text{max}}$ to the MH.

###  V.1 Boundary Conditions at the Universal Horizon and Spatial Infinity

To find the proper boundary conditions for ${\hat{\cal{Z}}}_{lm}(r)$ at $r\to
3/4$ and $r\to+\infty$, we first note that in terms of $x$ Eq.(3.32) takes the
form,

$\displaystyle\left[\frac{d^{2}}{dx^{2}}-V_{\text{eff}}(r)\right]\psi$
$\displaystyle=$ $\displaystyle 0.~{}~{}~{}~{}$ (5.2)

Note that (5.2) is at the same position as (3.27).

####  V.1.1 At Spatial Infinity

When $r\to+\infty$, (5.2) leads to $\psi={\bar{a}}_{+}+{\bar{a}}_{-}x$. For
${\hat{\cal{Z}}}_{lm}$ to be finite at $r\to+\infty$, we shall set
${\bar{a}}_{-}=0$ [cf. (3.33)]. Thus, we have

$\displaystyle\psi$ $\displaystyle=$
$\displaystyle\sum_{n=2}^{\infty}\frac{{\bar{a}}_{n}}{r^{n}},$ (5.3)

where ${\bar{a}}_{n}$ can be calculated through recursion relations similar to
(IV.1). Note that it is known a posteriori that the polynomial solution of
$\psi$ at $r\to+\infty$ starts from $n=2$, as already indicated in Eq.(5.3).
Combining with (3.33), we have found the polynomial solution of
${\hat{\cal{Z}}}_{lm}$ in the neighborhood of $r\to+\infty$.

####  V.1.2 At the Universal Horizon

On the other hand, at the UH $r=3/4$, by applying the formulas presented in
Appendix B, we can write (5.1) in the form,

$\displaystyle\left[{\bar{W}}_{g}(r)\frac{d^{2}}{d{\bar{y}}^{2}}+\left(-{\bar{V}}_{g}\right)\right]{\bar{\psi}}=0,$
(5.4)

where ${\bar{W}}_{g}\equiv\kappa_{1}/{\bar{p}}^{2}$, and

$\displaystyle{\bar{\psi}}$ $\displaystyle\equiv$
$\displaystyle-{\hat{\cal{Z}}}_{lm}\frac{\sqrt{\left(64r^{3}-16r^{2}-12r-9\right)\left(2\sqrt{2l^{2}+2l+9}+4r-3\right)}}{r^{2}}$
$\displaystyle\times\exp\left\\{\frac{1}{8}i\omega\left[3\sqrt{6}\tanh^{-1}\left(\frac{8r+3}{3\sqrt{8r^{2}+4r+\frac{3}{2}}}\right)-8\tanh^{-1}\left(\frac{20r+7}{3\sqrt{24r(2r+1)+9}}\right)\right]\right\\},$
$\displaystyle\frac{dr}{d{\bar{y}}}$ $\displaystyle\equiv$
$\displaystyle{\bar{p}}(r)=\frac{4r-3}{2\sqrt{2l^{2}+2l+9}+4r-3},$
$\displaystyle{\bar{V}}_{g}$ $\displaystyle\equiv$
$\displaystyle{-\left\\{\frac{{\bar{W}}_{g}}{4}\left[2{\bar{p}}\frac{d\left({\bar{p}}^{\prime}-{\bar{p}}\kappa_{2}/\kappa_{1}\right)}{dr}-\left({\bar{p}}^{\prime}-\frac{{\bar{p}}\kappa_{2}}{\kappa_{1}}\right)^{2}\right]+\kappa_{3}.\right\\}.}$
(5.5)

When $r\to(3/4)^{+}$, (5.4) leads to
${\bar{\psi}}={\bar{b}}_{+}e^{+{\bar{y}}}+{\bar{b}}_{-}e^{-{\bar{y}}}$. For
${\hat{\cal{Z}}}_{lm}$ to be finite at $r\to(3/4)^{+}$, we shall set
${\bar{b}}_{-}=0$ [cf. (3.33)]. Thus, we have

$\displaystyle{\bar{\psi}}$ $\displaystyle=$ $\displaystyle
e^{\bar{y}}\sum_{n=0}^{\infty}{\bar{b}}_{n}\left(r-\frac{3}{4}\right)^{n}.$
(5.6)

Combining with (3.33), we have found the polynomial solution of
${\hat{\cal{Z}}}_{lm}$ in the neighborhood of $r\to(3/4)^{+}$.

###  V.2 Smoothness Conditions Across the Metric Horizon

At the MH ($r=1$), it is clear that Eq.(5.1) becomes singular. By following
the same procedures of solving Eq.(4.13), we can also solve Eq.(5.1) in the
neighborhood of $r=1$. To do so, we shall first write
${\hat{\cal{Z}}}_{lm}(r)$ in the form,

$\displaystyle{\hat{\cal{Z}}}_{lm}(r)=(r-1)^{\bar{s}}\sum_{n=0}^{\infty}{\bar{d}}_{n}\left(r-1\right)^{n}.$
(5.7)

Just like how we obtain (4.18), the solution is determined to be of the form

$\displaystyle{\hat{\cal{Z}}}_{lm}(r)$ $\displaystyle=$
$\displaystyle(r-1)^{-i\omega}\left[\left(r-1\right)\sum_{n=0}^{\infty}{\bar{d}}_{n}^{-}\left(r-1\right)^{n}+{\cal{B}}\sum_{n=0}^{\infty}{\bar{d}}_{n}^{+}\left(r-1\right)^{n}\right],$
(5.8)

where ${\bar{d}}_{n}^{\pm}$ can be calculated through recursion relations
similar to (IV.3) and (4.17), and the (complex) constant ${\cal{B}}$ is acting
as ${\cal{A}}$ in (4.18).

###  V.3 Numerically Solving the aether Master Equation

To solve the aether master equation, let us first set

$\displaystyle{\hat{\cal{Z}}}_{lm}(r)$ $\displaystyle=$
$\displaystyle\text{Exp}\left(i\int^{r}\phi(r)dr\right),$ (5.9)

[It must not be confused with $\phi$ introduced here and the one appearing in
Eq.(2.24)] so that Eq.(5.1) takes the form,

$\alpha_{1}\phi^{\prime}(r)+\alpha_{2}\phi(r)^{2}+\alpha_{3}\phi(r)+\alpha_{4}=0,$
(5.10)

where $\alpha_{n}$’s are given in Eq.(Appendix A: The coefficients of
$\rho_{abc}$, $\kappa_{n}$ and $\alpha_{m}$) and must not be confused with the
PPN parameters mentioned in Sec. II. With the above given boundary conditions
at the universal horizon and spatial infinity, as well as the smoothness
conditions across the metric horizon, following the steps presented in Sec.
IV.D, we can solve Eq.(5.10) with the shooting method As what we did in the
last section, the new variable $z$ will be introduced in practice. In
addition, here in Table 3 we provide a summary to the usages of those main
equations that related to solving Eq.(5.10).

Table 3: Usages of equations in generating the suitable boundary conditions for solving Eq.(5.10). Type of boundary condition |  |  Polynomial |  |  Translated to $\phi(r)$ with
---|---|---|---|---
|  | expression |  |
Finity at $r\to+\infty$ |  | Eq.(5.3) |  | Eqs.(5.9), (3.33)
Finity at $r\to(3/4)^{+}$ |  | Eq.(5.6) |  | Eqs.(5.9), (V.1.2)
Smoothness across $r=1$ |  | Eq.(5.8) |  | Eq.(5.9)

However, different from what we have done in Sec. IV, in this section we would
like to address the problem of instability found in TZZW21  for large $l$. In
particular, in Table 4, in which $\Delta{\cal{B}}$ is defined in the same
manner as that of $\Delta{\cal{A}}$, we present the values of $\omega_{\ae}$
for $l=100$. In contrast to those given in Table 2, now the imaginary parts of
$\omega_{\ae}$’s are all positive, which indicates the instability of the
aether field. However, we must note the differences between the tolerance
levels of Table 2 indicated by $\Delta{\cal A}$ and the ones given in Table 4,
indicated by $\Delta\cal B$. In particular, the tolerance levels given by
$\Delta\cal B$ are much weaker than those given by $\Delta{\cal A}$.

Table 4: Values of $\omega_{\ae}$ by solving Eq.(5.10) for the aether field with $l=100$ and the units $c=G_{N}=r_{MH}=1$. $\omega_{\ae}$ |  | $\Delta{\cal B}$
---|---|---
$-1944.19+87.4996i$ |  |  $1.16\times 10^{-5}$
$-2913.80+155.310i$ |  |  $5.19\times 10^{-6}$
$-9886.49+267.356i$ |  |  $4.52\times 10^{-7}$
$-139279.4+101836.6i$ |  |  $1.13\times 10^{-8}$
$-264002.5+102057.4i$ |  |  $5.43\times 10^{-9}$

In Table 5, we also consider the cases $l=2,10$, from which it can be seen
that unstable modes in each case are also identified, again with much weaker
tolerance levels $\Delta{\cal B}\approx{\cal O}\left(10^{-10}\right)$.

Table 5: Values of $\omega_{\ae}$ by solving Eq.(5.10) for the aether field with the units $c=G_{N}=r_{MH}=1$ for $l=2,10$, respectively. $l$ | $\omega_{\ae}$ | $\Delta{\cal B}$
---|---|---
$2$ |  $-43880\pm 165.44i$ |  ${\cal O}\left(10^{-11}\right)~{}$
$10$ | $-46071\pm 291.48i$ | ${\cal O}\left(10^{-10}\right)$

##  VI Conclusions

In this paper, we have investigated QNMs of black holes in Einstein-aether
($\ae$-) theory, in which the gravitational sector of the theory is described
by three different species of gravitons with spins $0,1,2$, respectively JM01
. By properly choosing the four free dimensionless coupling constants,
$c_{i}$’s, the theory is self-consistent, such as free of ghosts and kinetic
instability Jacobson , and satisfies all the observations carried out so far
OMW18 ; GHBBCYY ; TZZW21 . However, to avoid the vacuum gravi-Čerenkov
radiation, such as cosmic rays, each of these three species must move with a
speed that is at least no less than the speed of light EMS05 .

Due to the existence of such superluminal motions, we expect that black holes
in $\ae$-theory will be dramatically different from those in GR. In
particular, the metric horizon that is used to define the boundary of a black
hole in GR is no longer a one-way membrane, and particles with superluminal
speeds can escape to the spatial infinity, even they are initially trapped
inside the metric horizon. As a result, the inner boundary of a black hole now
is defined by the universal horizon, a one-way membrane for particles with any
speeds, including the ones with infinitely large ones BS11 . This in turn
implies that at the metric horizon now both “in-going” and “out-going”
radiations are expected to exist, which will alternate the QNM spectra of of
black holes.

In addition, due to the presence of the aether field, the uniqueness theorem
of the Schwarzschild black hole does not hold any longer Jacobson , and
various black holes in $\ae$-theory exist Chao2020b , which can be still
formed from gravitational collapse of realistic matter Garfinkle2007 ,
including the formation of universal horizons BMSW18 .

The change of the location of a black hole inner boundary from the metric
horizon to the universal horizon also results in several technical challenges
in the studies of QNMs of black holes. First, the notion of “in-going” and
“out-going” radiations must be define properly, as now the universal horizon
is always inside the metric horizon, in which the Schwarzschild timelike
coordinate becomes spacelike. In addition, across the metric horizon, the
perturbation equations are singular, and proper smoothness conditions must be
imposed, too, before solving the corresponding boundary conditions imposed at
the universal horizon and the spatial infinity. In GR, the field equations of
perturbations are also singular at a metric horizon, but the problem is nicely
avoided, as now it is the inner boundary, and the smoothness conditions are
simply replaced by the inner boundary conditions.

To overcome the above problems, in this paper we have studied the QNMs of the
Schwarzschild black hole, which is also a solution of the Einstein-aether
theory Chao2020b , in which the presence of the aether field has no effects on
the spacetime, although the aether field is non-trivial [cf. Eqs.(II.2) and
(3.6)]. Despite the fact that the background is the same as that in GR, it is
expected that the corresponding QNMs are different, as now the inner boundary
is moved to the universal horizon, and both “in-going” and “out-going”
radiations exist at the metric horizon. In addition, the perturbation
equations are also different [cf. Eqs.(A. 5) and (3.23)].

In particular, in Section III we have constructed the two master equations
(3.27) and (3.32), respectively, for the metric and aether field
perturbations, while in Section IV we have first worked out step by step the
boundary conditions for the master equation of the metric perturbations at the
universal horizon, $r_{UH}=3r_{s}/4$, where $r_{s}=2m$ is the location of the
metric horizon, and at the spatial infinity. Then, we have systematically
worked out its smoothness conditions across the metric horizon. Afterwards,
combining the Chandrasekhar-Detweiler method Chandra1975  with the shooting
one, one side is from the universal horizon to the metric horizon, while the
other side is from the spatial infinity to the metric horizon [cf. Section
IV.D], we have found various modes $\omega$ for $l=2,3,4,10$, given explicitly
in Table 2. In the same table we have also presented the corresponding values
of the modes in GR, and found that the differences between the ones given in
$\ae$-theory and those given in GR can hardly be distinguished by current
observations of gravitational waves, although it is well within the
detectability of the detectors of the next generation.

In Section V, following closely what we have done for the metric perturbations
in Section IV, we have considered the aether perturbations and found various
modes $\omega_{\ae}$ [cf. Eq.(3.35)] for $l=2,10,100$, given, respectively in
Tables 4 and 5. In these two tables, in contrast to Table 2, we have
identified several modes of $\omega_{\ae}$’s, which all have a positive
imaginary part. This strongly indicates the instability of the aether field.
This is also consistent with what was found in TZZW21  for the Schwarzschild
black hole against large angular perturbations. However, due to our numerical
errors of the current paper, it must be noted that such a conclusion is not
exclusive. In particular, in the metric perturbations presented in Section IV,
we have not been able to find such an instability.

Therefore, it would be very important to confirm the above instability, for
example, by considering the even-parity perturbations and meanwhile extending
the current studies to other black holes found in Chao2020b . It would be also
very interesting to confirm the isospectrality of the theory Chao2022 .
However, due to the scope of this paper, we would like to address these
important questions in other occasions.

## Acknowledgments

We would like to express our gratitude to Prof. Shinji Tsujikawa for his
valuable comments and suggestions, and to Dr. Xiang Zhao for the early stage
of his collaboration of this project. This work is supported in part by the
National Key Research and Development Program of China Grant No.
2020YFC2201503, the National Natural Science Foundation of China under Grant
No. 12205254, 12275238, 11975203 and 11675143, the Zhejiang Provincial Natural
Science Foundation of China under Grant No. LR21A050001 and LY20A050002, and
the Fundamental Research Funds for the Provincial Universities of Zhejiang in
China under Grant No. RF-A2019015.

## Appendix A: The coefficients of $\rho_{abc}$, $\kappa_{n}$ and $\alpha_{m}$

With the gauge condition $C_{lm}=0$, the linearized Einstein-aether field
equations

$E_{\phi t}=E_{\phi r}=E_{\phi\theta}=\AE^{\phi}=0,$ (A. 1)

can be cast respectively in the forms,

$\displaystyle\rho_{101}{\ddot{J}}_{lm}+\rho_{102}{\dot{J}}_{lm}+\rho_{103}{\dot{J}}_{lm}^{\prime}+\rho_{104}a_{lm}+\rho_{105}{\ddot{a}}_{lm}+\rho_{106}{\dot{a}}_{lm}+\rho_{107}{\dot{a}}_{lm}^{\prime}+\rho_{108}a_{lm}^{\prime}+\rho_{109}a_{lm}^{\prime\prime}+\rho_{110}{\dot{G}}_{lm}=0,~{}~{}~{}~{}~{}~{}$
(A. 2)
$\displaystyle\rho_{201}{J}_{lm}+\rho_{202}{\dot{J}}_{lm}+\rho_{203}{\dot{J}}_{lm}^{\prime}+\rho_{204}a_{lm}+\rho_{205}{\ddot{a}}_{lm}+\rho_{206}{\dot{a}}_{lm}+\rho_{207}{\dot{a}}_{lm}^{\prime}+\rho_{208}a_{lm}^{\prime}+\rho_{209}a_{lm}^{\prime\prime}+\rho_{210}{G}_{lm}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\rho_{211}{G}_{lm}^{\prime}=0,~{}~{}~{}~{}~{}~{}~{}~{}$
(A. 3)
$\displaystyle\rho_{301}{J}_{lm}+\rho_{302}{\dot{J}}_{lm}+\rho_{303}{J}_{lm}^{\prime}+\rho_{304}G_{lm}+\rho_{305}{\dot{G}}_{lm}+\rho_{306}{\dot{G}}_{lm}^{\prime}+\rho_{307}G_{lm}^{\prime}+\rho_{308}G_{lm}^{\prime\prime}=0,~{}~{}~{}~{}~{}~{}~{}~{}$
(A. 4)
$\displaystyle\rho_{401}a_{lm}+\rho_{402}{\ddot{a}}_{lm}+\rho_{403}{\dot{a}}_{lm}+\rho_{404}{\dot{a}}_{lm}^{\prime}+\rho_{405}a_{lm}^{\prime}+\rho_{406}a_{lm}^{\prime\prime}=0,~{}~{}~{}~{}~{}~{}~{}~{}$
(A. 5)

where

$\displaystyle\rho_{101}$ $\displaystyle=$ $\displaystyle-\frac{1}{2},$
$\displaystyle\rho_{102}$ $\displaystyle=$ $\displaystyle\frac{1-r}{r^{2}},$
$\displaystyle\rho_{103}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{1}{r}-1\right),$
$\displaystyle\rho_{104}$ $\displaystyle=$
$\displaystyle\left\\{256r^{8}\left(16r^{2}+8r+3\right)\sqrt{256(r-1)r^{3}+27}\left[\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right]^{5}\right\\}^{-1}$
$\displaystyle\times\Big{\\{}c_{1}\Big{[}98304\sqrt{3}r^{8}-196608\sqrt{3}r^{7}+98304\sqrt{3}r^{6}+31104\sqrt{3}r^{4}-31104\sqrt{3}r^{3}+6912\sqrt{256(r-1)r^{3}+27}r^{3}$
$\displaystyle-729\sqrt{256(r-1)r^{3}+27}-8192\sqrt{256(r-1)r^{3}+27}r^{8}+16384\sqrt{256(r-1)r^{3}+27}r^{7}$
$\displaystyle-8192\sqrt{256(r-1)r^{3}+27}r^{6}-6912\sqrt{256(r-1)r^{3}+27}r^{4}+2187\sqrt{3}\Big{]}$
$\displaystyle\times\left[-256r^{4}+256r^{3}+3\sqrt{768(r-1)r^{3}+81}-27\right]$
$\displaystyle\times\left[1024\left(l^{2}+l-2\right)r^{5}+128\left(3l^{2}+3l-4\right)r^{4}+2048l(l+1)r^{6}+384r^{3}+3024r^{2}+1512r+567\right]\Big{\\}},$
$\displaystyle\rho_{105}$ $\displaystyle=$
$\displaystyle-\frac{c_{1}\left[3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right]\left[-256r^{4}+256r^{3}+3\sqrt{768(r-1)r^{3}+81}-27\right]}{4096(r-1)^{2}r^{4}},$
$\displaystyle\rho_{106}$ $\displaystyle=$
$\displaystyle\frac{3c_{1}(4r-3)\left[\sqrt{768(r-1)r^{3}+81}-9\right]\left[-256r^{4}+256r^{3}+3\sqrt{768(r-1)r^{3}+81}-27\right]}{8192(r-1)^{2}r^{6}\sqrt{256(r-1)r^{3}+27}},$
$\displaystyle\rho_{107}$ $\displaystyle=$
$\displaystyle\frac{c_{1}\left[256r^{4}-256r^{3}-3\sqrt{768(r-1)r^{3}+81}+27\right]^{2}}{2048(r-1)r^{5}\left[3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right]},$
$\displaystyle\rho_{108}$ $\displaystyle=$
$\displaystyle\Big{\\{}512r^{7}\sqrt{256(r-1)r^{3}+27}\left[\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right]^{4}\Big{\\}}^{-1}$
$\displaystyle\times
c_{1}\left(64r^{3}-27\right)\left[256r^{4}-256r^{3}-3\sqrt{768(r-1)r^{3}+81}+27\right]^{2}$
$\displaystyle\times\left[-128r^{4}+128r^{3}+3\sqrt{768(r-1)r^{3}+81}-27\right],$
$\displaystyle\rho_{109}$ $\displaystyle=$
$\displaystyle-\frac{c_{1}\left[-256r^{4}+256r^{3}+3\sqrt{768(r-1)r^{3}+81}-27\right]^{3}}{4096r^{6}\left[3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right]^{3}},$
$\displaystyle\rho_{110}$ $\displaystyle=$
$\displaystyle-\frac{l^{2}+l-2}{4r^{2}},$ (A. 6) $\displaystyle\rho_{201}$
$\displaystyle=$ $\displaystyle-\frac{l^{2}+l-2}{2r^{2}},$
$\displaystyle\rho_{202}$ $\displaystyle=$ $\displaystyle\frac{1}{r},$
$\displaystyle\rho_{203}$ $\displaystyle=$ $\displaystyle\frac{1}{2},$
$\displaystyle\rho_{204}$ $\displaystyle=$
$\displaystyle\Big{\\{}256(r-1)r^{7}\left(16r^{2}+8r+3\right)\sqrt{256(r-1)r^{3}+27}\left[\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right]^{3}\Big{\\}}^{-1}$
$\displaystyle\times
c_{1}\Big{[}98304\sqrt{3}r^{8}-196608\sqrt{3}r^{7}+98304\sqrt{3}r^{6}+31104\sqrt{3}r^{4}-31104\sqrt{3}r^{3}+6912\sqrt{256(r-1)r^{3}+27}r^{3}$
$\displaystyle-729\sqrt{256(r-1)r^{3}+27}-8192\sqrt{256(r-1)r^{3}+27}r^{8}+16384\sqrt{256(r-1)r^{3}+27}r^{7}$
$\displaystyle-8192\sqrt{256(r-1)r^{3}+27}r^{6}-6912\sqrt{256(r-1)r^{3}+27}r^{4}+2187\sqrt{3}\Big{]}$
$\displaystyle\times\left[1024\left(l^{2}+l-2\right)r^{5}+128\left(3l^{2}+3l-4\right)r^{4}+2048l(l+1)r^{6}+384r^{3}+3024r^{2}+1512r+567\right],$
$\displaystyle\rho_{205}$ $\displaystyle=$
$\displaystyle-\frac{c_{1}\left[3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right]^{3}}{4096(r-1)^{3}r^{3}},$
$\displaystyle\rho_{206}$ $\displaystyle=$
$\displaystyle\frac{3c_{1}(4r-3)\left[\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right]^{2}\left[\sqrt{768(r-1)r^{3}+81}-9\right]}{8192(r-1)^{3}r^{5}\sqrt{256(r-1)r^{3}+27}},$
$\displaystyle\rho_{207}$ $\displaystyle=$
$\displaystyle\frac{c_{1}\left[3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right]\left[-256r^{4}+256r^{3}+3\sqrt{768(r-1)r^{3}+81}-27\right]}{2048(r-1)^{2}r^{4}},$
$\displaystyle\rho_{208}$ $\displaystyle=$
$\displaystyle\frac{c_{1}\left(64r^{3}-27\right)\left(256r^{4}-256r^{3}-3\sqrt{768r^{4}-768r^{3}+81}+27\right)}{1024(r-1)r^{6}\sqrt{256r^{4}-256r^{3}+27}},$
$\displaystyle\rho_{209}$ $\displaystyle=$
$\displaystyle-\frac{c_{1}\left[256r^{4}-256r^{3}-3\sqrt{768(r-1)r^{3}+81}+27\right]^{2}}{4096(r-1)r^{5}\left[3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right]},$
$\displaystyle\rho_{210}$ $\displaystyle=$
$\displaystyle\frac{l^{2}+l-2}{2r^{3}},$ $\displaystyle\rho_{211}$
$\displaystyle=$ $\displaystyle-\frac{l^{2}+l-2}{4r^{2}},$ (A. 7)
$\displaystyle\rho_{301}$ $\displaystyle=$ $\displaystyle\frac{1}{2r^{2}},$
$\displaystyle\rho_{302}$ $\displaystyle=$ $\displaystyle\frac{1}{2},$
$\displaystyle\rho_{303}$ $\displaystyle=$ $\displaystyle\frac{r-1}{2r},$
$\displaystyle\rho_{304}$ $\displaystyle=$ $\displaystyle\frac{r-2}{2r^{3}},$
$\displaystyle\rho_{305}$ $\displaystyle=$ $\displaystyle-\frac{1}{2r},$
$\displaystyle\rho_{306}$ $\displaystyle=$ $\displaystyle\frac{1}{2},$
$\displaystyle\rho_{307}$ $\displaystyle=$ $\displaystyle\frac{3-2r}{4r^{2}},$
$\displaystyle\rho_{308}=$ $\displaystyle=$ $\displaystyle\frac{r-1}{4r},$ (A.
8) $\displaystyle\rho_{401}$ $\displaystyle=$
$\displaystyle\Big{\\{}16r^{8}\left(16r^{2}+8r+3\right)\sqrt{256(r-1)r^{3}+27}\left[\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right]^{4}\Big{\\}}^{-1}$
$\displaystyle\times\Big{\\{}\Big{[}-98304\sqrt{3}r^{8}+196608\sqrt{3}r^{7}-98304\sqrt{3}r^{6}-31104\sqrt{3}r^{4}+31104\sqrt{3}r^{3}-6912\sqrt{256(r-1)r^{3}+27}r^{3}$
$\displaystyle+729\sqrt{256(r-1)r^{3}+27}+8192\sqrt{256(r-1)r^{3}+27}r^{8}-16384\sqrt{256(r-1)r^{3}+27}r^{7}$
$\displaystyle+8192\sqrt{256(r-1)r^{3}+27}r^{6}+6912\sqrt{256(r-1)r^{3}+27}r^{4}-2187\sqrt{3}\Big{]}$
$\displaystyle\times\left[1024\left(l^{2}+l-2\right)r^{5}+128\left(3l^{2}+3l-4\right)r^{4}+2048l(l+1)r^{6}+384r^{3}+3024r^{2}+1512r+567\right]\Big{\\}},$
$\displaystyle\rho_{402}$ $\displaystyle=$
$\displaystyle-\frac{\left(\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right)^{2}}{256(r-1)^{2}r^{4}},$
$\displaystyle\rho_{403}$ $\displaystyle=$
$\displaystyle-\frac{3(4r-3)\left[\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right]\left[\sqrt{768(r-1)r^{3}+81}-9\right]}{512(r-1)^{2}r^{6}\sqrt{256(r-1)r^{3}+27}},$
$\displaystyle\rho_{404}$ $\displaystyle=$
$\displaystyle\frac{-256r^{4}+256r^{3}+3\left(\sqrt{768r^{4}-768r^{3}+81}-9\right)}{128(r-1)r^{5}},$
$\displaystyle\rho_{405}$ $\displaystyle=$
$\displaystyle-\frac{\left(64r^{3}-27\right)\left[128r^{4}-128r^{3}-3\sqrt{768(r-1)r^{3}+81}+27\right]\left[256r^{4}-256r^{3}-3\sqrt{768(r-1)r^{3}+81}+27\right]}{32r^{7}\sqrt{256(r-1)r^{3}+27}\left(\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right)^{3}},$
$\displaystyle\rho_{406}$ $\displaystyle=$
$\displaystyle-\frac{27}{256r^{6}}+\frac{1}{r^{3}}-\frac{1}{r^{2}}.$ (A. 9)

It is interesting to note that the coefficients $\rho_{abc}$ are functions of
$r$, $c_{1}$ and $l$ only, and some of them contain no $c_{1}$.

The coefficients $\kappa_{n}$ appearing in Eq.(5.1) are given by

$\displaystyle\kappa_{0}$ $\displaystyle\equiv$
$\displaystyle\sqrt{256(r-1)r^{3}+27},$ $\displaystyle\kappa_{1}$
$\displaystyle\equiv$ $\displaystyle\frac{(r-1)^{2}}{r^{2}},$
$\displaystyle\kappa_{2}$ $\displaystyle\equiv$
$\displaystyle\Big{[}4194304r^{14}-9437184r^{13}+5505024r^{12}-65536r^{11}+2916\left(\sqrt{768(r-1)r^{3}+81}-9\right)r^{3}$
(A. 10)
$\displaystyle-49152\left(3\sqrt{768(r-1)r^{3}+81}-53\right)r^{10}+36864\left(5\sqrt{768(r-1)r^{3}+81}-91\right)r^{9}$
$\displaystyle-9216\left(\sqrt{768(r-1)r^{3}+81}-15\right)r^{8}-6912\left(\sqrt{768(r-1)r^{3}+81}-15\right)r^{7}$
$\displaystyle-41472\left(\sqrt{768(r-1)r^{3}+81}-12\right)r^{6}+5184\left(\sqrt{768(r-1)r^{3}+81}-9\right)r^{5}$
$\displaystyle+3888\left(\sqrt{768(r-1)r^{3}+81}-9\right)r^{4}\Big{]}^{-1}$
$\displaystyle\times\Big{\\{}\Big{[}4194304r^{12}-13631488r^{11}+14942208r^{10}-7340032r^{9}+15552\left(7\sqrt{768(r-1)r^{3}+81}-95\right)r^{3}$
$\displaystyle+7776\left(\sqrt{768(r-1)r^{3}+81}-9\right)r-11664\left(\sqrt{768(r-1)r^{3}+81}-9\right)$
$\displaystyle-16384\left(9\sqrt{768(r-1)r^{3}+81}-487\right)r^{8}+12288\left(27\sqrt{768(r-1)r^{3}+81}-917\right)r^{7}$
$\displaystyle-9216\left(21\sqrt{768(r-1)r^{3}+81}-571\right)r^{6}+64512\left(\sqrt{768(r-1)r^{3}+81}-15\right)r^{5}$
$\displaystyle-6912\left(23\sqrt{768(r-1)r^{3}+81}-327\right)r^{4}-1296\left(\sqrt{768(r-1)r^{3}+81}-9\right)r^{2}\Big{]}$
$\displaystyle+\Big{[}1944i\left(\sqrt{768(r-1)r^{3}+81}-9\right)r^{3}+98304i\left(\sqrt{768(r-1)r^{3}+81}-27\right)r^{10}$
$\displaystyle-221184i\left(\sqrt{768(r-1)r^{3}+81}-27\right)r^{9}+129024i\left(\sqrt{768(r-1)r^{3}+81}-27\right)r^{8}$
$\displaystyle-1536i\left(\sqrt{768(r-1)r^{3}+81}-27\right)r^{7}+4608i\left(11\sqrt{768(r-1)r^{3}+81}-135\right)r^{6}$
$\displaystyle-3456i\left(19\sqrt{768(r-1)r^{3}+81}-243\right)r^{5}+2592i\left(\sqrt{768(r-1)r^{3}+81}-9\right)r^{4}$
$\displaystyle+5832i\left(\sqrt{768(r-1)r^{3}+81}-9\right)r^{2}\Big{]}\omega_{\ae}\Big{\\}},$
$\displaystyle\kappa_{3}$ $\displaystyle\equiv$
$\displaystyle\Big{[}-268435456\left(12\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right)r^{18}+671088640\left(12\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right)r^{17}$
(A. 11)
$\displaystyle-452984832\left(12\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right)r^{16}+16777216\left(12\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right)r^{15}$
$\displaystyle-1048576\left(1488\sqrt{3}-259\sqrt{256(r-1)r^{3}+27}\right)r^{14}+7864320\left(336\sqrt{3}-55\sqrt{256(r-1)r^{3}+27}\right)r^{13}$
$\displaystyle-15925248\left(16\sqrt{3}-3\sqrt{256(r-1)r^{3}+27}\right)r^{12}-10616832\left(16\sqrt{3}-3\sqrt{256(r-1)r^{3}+27}\right)r^{11}$
$\displaystyle-3981312\left(109\sqrt{3}-24\sqrt{256(r-1)r^{3}+27}\right)r^{10}+5971968\left(15\sqrt{3}-4\sqrt{256(r-1)r^{3}+27}\right)r^{9}$
$\displaystyle+3732480\left(15\sqrt{3}-4\sqrt{256(r-1)r^{3}+27}\right)r^{8}+2239488\left(15\sqrt{3}-4\sqrt{256(r-1)r^{3}+27}\right)r^{7}$
$\displaystyle-2519424\left(3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right)r^{6}-1259712\left(3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right)r^{5}$
$\displaystyle-472392\left(3\sqrt{3}-\sqrt{256(r-1)r^{3}+27}\right)r^{4}\Big{]}^{-1}$
$\displaystyle\times\Big{\\{}\Big{[}-268435456\left(12\sqrt{3}-\kappa_{0}\right)r^{15}+1006632960\left(12\sqrt{3}-\kappa_{0}\right)r^{14}-1291845632\left(12\sqrt{3}-\kappa_{0}\right)r^{13}$
$\displaystyle+866123776\left(12\sqrt{3}-\kappa_{0}\right)r^{12}-9437184\left(1488\sqrt{3}-139\kappa_{0}\right)r^{11}+1441792\left(13524\sqrt{3}-1397\kappa_{0}\right)r^{10}$
$\displaystyle-786432\left(12489\sqrt{3}-1412\kappa_{0}\right)r^{9}+2654208\left(671\sqrt{3}-118\kappa_{0}\right)r^{8}-1327104\left(3161\sqrt{3}-584\kappa_{0}\right)r^{7}$
$\displaystyle+497664\left(6763\sqrt{3}-1282\kappa_{0}\right)r^{6}-1492992\left(78\sqrt{3}-19\kappa_{0}\right)r^{5}+186624\left(831\sqrt{3}-221\kappa_{0}\right)r^{4}$
$\displaystyle-279936\left(1527\sqrt{3}-397\kappa_{0}\right)r^{3}+6613488\left(3\sqrt{3}-\kappa_{0}\right)r^{2}+4408992\left(3\sqrt{3}-\kappa_{0}\right)r\Big{]}$
$\displaystyle+\left[28311552\left(12\sqrt{3}-\kappa_{0}\right)r^{14}-42467328\left(12\sqrt{3}-\kappa_{0}\right)r^{13}+5308416\left(12\sqrt{3}-\kappa_{0}\right)r^{12}\right.$
$\displaystyle\left.+3538944\left(12\sqrt{3}-\kappa_{0}\right)r^{11}+1327104\left(129\sqrt{3}-22\kappa_{0}\right)r^{10}-5971968\left(9\sqrt{3}-2\kappa_{0}\right)r^{9}\right.$
$\displaystyle\left.-3732480\left(9\sqrt{3}-2\kappa_{0}\right)r^{8}-2239488\left(9\sqrt{3}-2\kappa_{0}\right)r^{7}+2519424\left(3\sqrt{3}-\kappa_{0}\right)r^{6}+1259712\left(3\sqrt{3}-\kappa_{0}\right)r^{5}\right.$
$\displaystyle\left.+472392\left(3\sqrt{3}-\kappa_{0}\right)r^{4}\right]\omega_{\ae}^{2}+6613488\left(3\sqrt{3}-\kappa_{0}\right)$
$\displaystyle+l^{2}\Big{[}268435456\left(12\sqrt{3}-\kappa_{0}\right)r^{16}-939524096\left(12\sqrt{3}-\kappa_{0}\right)r^{15}+1124073472\left(12\sqrt{3}-\kappa_{0}\right)r^{14}$
$\displaystyle-469762048\left(12\sqrt{3}-\kappa_{0}\right)r^{13}+4194304\left(339\sqrt{3}-62\kappa_{0}\right)r^{12}-2097152\left(1599\sqrt{3}-302\kappa_{0}\right)r^{11}$
$\displaystyle-+3932160\left(591\sqrt{3}-110\kappa_{0}\right)r^{10}-7077888\left(9\sqrt{3}-2\kappa_{0}\right)r^{9}+884736\left(153\sqrt{3}-43\kappa_{0}\right)r^{8}$
$\displaystyle-1327104\left(225\sqrt{3}-59\kappa_{0}\right)r^{7}+4478976\left(3\sqrt{3}-\kappa_{0}\right)r^{6}+2985984\left(3\sqrt{3}-\kappa_{0}\right)r^{5}+4478976\left(3\sqrt{3}-\kappa_{0}\right)r^{4}\Big{]}$
$\displaystyle+l\Big{[}268435456\left(12\sqrt{3}-\kappa_{0}\right)r^{16}-939524096\left(12\sqrt{3}-\kappa_{0}\right)r^{15}+1124073472\left(12\sqrt{3}-\kappa_{0}\right)r^{14}$
$\displaystyle-469762048\left(12\sqrt{3}-\kappa_{0}\right)r^{13}+4194304\left(339\sqrt{3}-62\kappa_{0}\right)r^{12}-2097152\left(1599\sqrt{3}-302\kappa_{0}\right)r^{11}$
$\displaystyle+3932160\left(591\sqrt{3}-110\kappa_{0}\right)r^{10}-7077888\left(9\sqrt{3}-2\kappa_{0}\right)r^{9}+884736\left(153\sqrt{3}-43\kappa_{0}\right)r^{8}$
$\displaystyle-1327104\left(225\sqrt{3}-59\kappa_{0}\right)r^{7}+4478976\left(3\sqrt{3}-\kappa_{0}\right)r^{6}+2985984\left(3\sqrt{3}-\kappa_{0}\right)r^{5}+4478976\left(3\sqrt{3}-\kappa_{0}\right)r^{4}\Big{]}$
$\displaystyle+\Big{[}-1610612736i\sqrt{3}r^{17}+5234491392i\sqrt{3}r^{16}-5737807872i\sqrt{3}r^{15}+2139095040i\sqrt{3}r^{14}$
$\displaystyle-12582912i\left(149\sqrt{3}-18\kappa_{0}\right)r^{13}+9437184i\left(445\sqrt{3}-54\kappa_{0}\right)r^{12}-7077888i\left(347\sqrt{3}-42\kappa_{0}\right)r^{11}$
$\displaystyle+1769472i\left(15\sqrt{3}-2\kappa_{0}\right)r^{10}-5308416i\left(111\sqrt{3}-22\kappa_{0}\right)r^{9}+3981312i\left(195\sqrt{3}-38\kappa_{0}\right)r^{8}$
$\displaystyle-2985984i\left(9\sqrt{3}-2\kappa_{0}\right)r^{7}-2239488i\left(9\sqrt{3}-2\kappa_{0}\right)r^{6}-3359232i\left(27\sqrt{3}-7\kappa_{0}\right)r^{5}$
$\displaystyle+2519424i\left(3\sqrt{3}-\kappa_{0}\right)r^{4}+1889568i\left(3\sqrt{3}-\kappa_{0}\right)r^{3}+1417176i\left(3\sqrt{3}-\kappa_{0}\right)r^{2}\Big{]}\omega_{\ae}\Big{\\}}.$

The coefficients $\alpha_{n}$ appearing in Eq.(5.10) are given by

$\displaystyle\alpha_{1}$ $\displaystyle=$
$\displaystyle\frac{i(r-1)^{2}}{r^{2}},$ $\displaystyle\alpha_{2}$
$\displaystyle=$ $\displaystyle-\frac{(r-1)^{2}}{r^{2}},$
$\displaystyle\alpha_{3}$ $\displaystyle=$
$\displaystyle-\frac{6(r-1)\sqrt{768r^{4}-768r^{3}+81}\omega_{\text{\ae}}}{256r^{5}-256r^{4}+27r}$
$\displaystyle+\frac{16i(r-1)^{2}\left(16r^{2}+12r+9\right)}{r^{3}(4r-3)\left(16r^{2}+8r+3\right)\sqrt{256(r-1)r^{3}+27}\left[\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right]^{3}}$
$\displaystyle\times\Big{[}16384r^{8}-32768r^{7}+16384r^{6}+8640r^{4}+576\sqrt{768(r-1)r^{3}+81}r^{3}-8640r^{3}$
$\displaystyle~{}~{}~{}-81\sqrt{768(r-1)r^{3}+81}-576\sqrt{768(r-1)r^{3}+81}r^{4}+729\Big{]},$
$\displaystyle\alpha_{4}$ $\displaystyle=$
$\displaystyle-\frac{27\omega_{\text{\ae}}^{2}}{256r^{4}-256r^{3}+27}$ (A. 12)
$\displaystyle-\frac{12i\left(128r^{4}-128r^{3}+27\right)\omega_{\text{\ae}}}{r^{2}(4r-3)\left(16r^{2}+8r+3\right)\sqrt{256(r-1)r^{3}+27}\left(\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right)^{3}}$
$\displaystyle\times\left[-1728r^{4}-64\sqrt{3}\kappa_{0}r^{3}+1728r^{3}+27\sqrt{3}\kappa_{0}+64\sqrt{3}\kappa_{0}r^{4}-243\right]$
$\displaystyle+\frac{16(r-1)^{2}\left(1024\left(l^{2}+l-2\right)r^{5}+128\left(3l^{2}+3l-4\right)r^{4}+2048l(l+1)r^{6}+384r^{3}+3024r^{2}+1512r+567\right)}{(3-4r)^{2}r^{4}\left(16r^{2}+8r+3\right)^{2}\sqrt{256(r-1)r^{3}+27}\left(\sqrt{256(r-1)r^{3}+27}-3\sqrt{3}\right)^{4}}$
$\displaystyle\times\Big{(}98304\sqrt{3}r^{8}-196608\sqrt{3}r^{7}+98304\sqrt{3}r^{6}+31104\sqrt{3}r^{4}-31104\sqrt{3}r^{3}+6912\kappa_{0}r^{3}-729\kappa_{0}-8192\kappa_{0}r^{8}$
$\displaystyle~{}~{}~{}+16384\kappa_{0}r^{7}-8192\kappa_{0}r^{6}-6912\kappa_{0}r^{4}+2187\sqrt{3}\Big{)}.$

## Appendix B: The formulas used in modifying a master equation

Suppose we have an eqaution of the form

$\displaystyle\left[\eta_{1}(r)\frac{d^{2}}{dr^{2}}+\eta_{2}(r)\frac{d}{dr}+\eta_{3}(r)\right]J(r)=0.$
(B. 1)

Let us introduce a new function $\Psi$ and a new variable $x$ by

$\displaystyle J=q(r)\Psi,\quad\frac{dr}{dx}=p(r).$ (B. 2)

Then, choosing $q(r)$ as

$\displaystyle q$ $\displaystyle=$
$\displaystyle\sqrt{p(r)}\exp{\left(-\frac{1}{2}\int\frac{\eta_{2}}{\eta_{1}}dr\right)},$
(B. 3)

we find that Eq.(B. 1) can be written in a form of the Schrödinger-like
differential equation

$\displaystyle\left[\frac{\eta_{1}}{p^{2}}\frac{d^{2}}{dx^{2}}-V(r)\right]\Psi=0,$
(B. 4)

where

$\displaystyle V(r)$ $\displaystyle\equiv$
$\displaystyle-\Bigg{\\{}\frac{\eta_{1}}{4p^{2}}\bigg{[}2p\frac{d\left(p^{\prime}-p\eta_{2}/\eta_{1}\right)}{dr}$
(B. 5)
$\displaystyle~{}~{}~{}~{}~{}~{}-\left(p^{\prime}-p\frac{\eta_{2}}{\eta_{1}}\right)^{2}\bigg{]}+\eta_{3}\Bigg{\\}},$

and a prime denotes the derivative with respect to $r$. Note that another
useful form of Eq.(B. 4) is given by

$\displaystyle\left[\frac{\eta_{1}}{p^{2}}\left(p^{2}\frac{d^{2}}{dr^{2}}+pp^{\prime}\frac{d}{dr}\right)-V(r)\right]\Psi=0,$
(B. 6)

where Eq.(B. 2) has been used. In addition, a common choice of $p(r)$ is
$p(r)=\sqrt{\eta_{1}}$.

## References

  * (1) B.P. Abbott, et al., [LIGO/Virgo Scientific Collaborations], Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016).
  * (2) B.P. Abbott, et al., [LIGO/Virgo Collaborations], GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X9, 031040 (2019).
  * (3) B.P. Abbott, et al., [LIGO/Virgo Collaborations], Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo, SoftwareX, Volume 13, 100658 (2021).
  * (4) B.P. Abbott, et al., [LIGO/Virgo Collaborations], GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\bigodot}$, ApJL 892 L3 (2020).
  * (5) B.P. Abbott, et al., [LIGO/Virgo/KAGRA Collaborations], GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, arXiv:2111.03606v1 [gr-qc].
  * (6) C. J. Moore, R. H. Cole and C. P. L. Berry, Gravitational-wave sensitivity curves, Class. Quantum. Grav. 32, 015014 (2015).
  * (7) Y. Gong, J. Luo and B. Wang, “Concepts and status of Chinese space gravitational wave detection projects,” Nature Astron. 5, no.9, 881-889 (2021).
  * (8) E. Berti, V. Cardoso and A. O. Starinets, Quasinormal modes of black holes and black branes, Class. Quantum. Grav. 26, 163001 (2009).
  * (9) E. Berti, K. Yagi, H. Yang, N. Yunes, Extreme gravity tests with gravitational waves from compact binary coalescences: (II) ringdown, Gen. Relativ. Grav. 50, 49 (2018).
  * (10) E. Berti, K. Yagi, H. Yang, N. Yunes, Extreme gravity tests with gravitational waves from compact binary coalescences: (I) inspiral-merger, Gen. Relativ. Grav. 50, 46 (2018).
  * (11) CE, https://cosmicexplorer.org/.
  * (12) ET Steering Committee Editorial Team, ET design report update 2020, ET-0007A- 20 (2020); https://www.et-gw.eu/.
  * (13) https://www.lisamission.org
  * (14) S. Liu, Y. Hu, et al., Science with the TianQin observatory: Preliminary results on stellar-mass binary black holes, Phys. Rev. D101, 103027 (2020); C.-F. Shi, et al., Science with the TianQin observatory: Preliminary results on testing the no-hair theorem with ringdown signals, Phys. Rev. D100, 044036 (2019).
  * (15) W.-H. Ruan, Z.-K. Guo, R.-G. Cai, Y.-Z. Zhang, Taiji Program: Gravitational-Wave Sources, Int. J. Mod. Phys. A 35, No. 17, 2050075 (2020).
  * (16) S. Kawamura, et al., Current status of space gravitational wave antenna DECIGO and B-DECIGO, arXiv:2006.13545.
  * (17) S. Chandrasekhar, the mathematical theory of black holes, Oxford classic texts in the physical sciences (Oxford Press, Oxford, 1992).
  * (18) S. Iyer, Black-hole normal modes: A WKB approach. II. Schwarzschild black holes, Phys. Rev. D35, 12 (1987).
  * (19) S. Detweiler, BLACK HOLES AND GRAVITATIONAL WAVES. III. THE RESONANT FREQUENCIES OF ROTATING HOLES, Astrophys. J. 239, 292-295 (1980).
  * (20) E. Seidel and S. Iyer, Black-hole normal modes: A WKB approach. IV. Kerr black holes, Phys. Rev. D41, 2 (1990).
  * (21) B. F. Schutz and C. M. Will, BLACK HOLE NORMAL MODES: A SEMIANALYTIC APPROACH, Astrophys. J. 291, L33-L36 (1985).
  * (22) S. Iyer and C. M. Will, Black-hole normal modes: A WKB approach. I. Foundations and application of a higher-order WKB analysis of potential-barrier scattering, Phys. Rev. D35, 12 (1987).
  * (23) R. A. Konoplya, Quasinormal behavior of the D-dimensional Schwarzschild black hole and the higher order WKB approach, Phys. Rev. D68, 024018 (2003).
  * (24) J. Matyjasek and M. Opala, Quasinormal modes of black holes: The improved semianalytic approach, Phys. Rev. D96, 024011 (2017).
  * (25) X. Li and S.-P. Zhao, Quasinormal modes of a scalar and an electromagnetic field in Finslerian-Schwarzschild spacetime, Phys. Rev. D101, 124012 (2020).
  * (26) E. W. Leaver, An analytic representation for the quasi-normal modes of Kerr black holes, Proc. R. Soc. Lond. A. 402, 285-298 (1985).
  * (27) S. Chandrasekhar, F. R. S., and S. Detweiler, The quasi-normal modes of the Schwarzschild black hole, Proc. R. Soc. Lond. A. 344, 411-452 (1975).
  * (28) D. D. Doneva, S. S. Yazadjiev, K. D. Kokkotas, and I. Zh. Stefanov, Quasinormal modes, bifurcations, and nonuniqueness of charged scalar-tensor black holes, Phys. Rev. D82, 064030 (2010).
  * (29) K. Lin and W.-L. Qian, A matrix method for quasinormal modes: Schwarzschild black holes in asymptotically flat and (anti-) de Sitter spacetimes, Class. Quantum Grav. 34, 095004 (2017).
  * (30) R. A. Konoplya and A. Zhidenko, Quasinormal modes of black holes: From astrophysics to string theory, Rev. Mod. Phys. 83, 793 (2011).
  * (31) C. Gundlach, R. H. Price and J. Pullin, Late-time behavior of stellar collapse and explosions. I. Linearized perturbations, Phys. Rev. D49, 883 (1994).
  * (32) B. Wang, C.-Y. Lin and C. Molina, Quasinormal behavior of massless scalar field perturbation in Reissner-Nordström anti-de Sitter spacetimes, Phys. Rev. D70, 064025 (2004).
  * (33) X. Li and S.-P. Zhao, Quasinormal modes of a scalar and an electromagnetic field in Finslerian-Schwarzschild spacetime, Phys. Rev. D101, 124012 (2020).
  * (34) O. J. Tattersall and P. G. Ferreira, Forecasts for low spin black hole spectroscopy in Horndeski gravity, Phys. Rev. D99, 104082 (2019).
  * (35) H. Yang, D. A. Nichols, F. Zhang, et al., Quasinormal-mode spectrum of Kerr black holes and its geometric interpretation, Phys. Rev. D86, 104006 (2012).
  * (36) T. Jacobson and D. Mattingly, Gravity with a dynamical preferred frame, Phys. Rev. D64, 024028 (2001).
  * (37) T. Jacobson, Einstein-$\ae$ther gravity: a status report, Proc. Sci., QG-PH (2007) 020 [arXiv:0801.1547v2].
  * (38) J. Oost, S. Mukohyama and A. Wang, Constraints on $\ae$-theory after GW170817, Phys. Rev. D97, 124023 (2018).
  * (39) O. Sarbach, E. Barausse, and J.A Preciado-López, Well-posed Cauchy formulation for Einstein-$\ae$ther theory, Class. Quantum Grav. 36 (2019) 165007.
  * (40) D. Garfinkle and T. Jacobson, A Positive-Energy Theorem for Einstein-Aether and Hořava Gravity, Phys. Rev. Lett. 107, 191102 (2011).
  * (41) B. Z. Foster, Radiation Damping in Einstein-Aether Theory, arXiv:gr-qc/0602004v5.
  * (42) B. Z. Foster, Phys. Rev. D76, 084033 (2007).
  * (43) K. Yagi, D. Blas, N. Yunes and E. Barausse, “Strong Binary Pulsar Constraints on Lorentz Violation in Gravity,” Phys. Rev. Lett. 112, 161101 (2014).
  * (44) K. Yagi, D. Blas, E. Barausse, and N. Yunes, Phys. Rev. D89, 084067 (2014).
  * (45) D. Hansen, N. Yunes, and K. Yagi, Phys. Rev. D91, 082003 (2015).
  * (46) Y.-G. Gong, S.-Q. Hou, D.-C. Liang, E. Papantonopoulos, Phys. Rev. D97, 084040 (2018).
  * (47) K. Lin, X. Zhao, C. Zhang, K. Lin, T. Liu, B. Wang, S.-J. Zhang, X. Zhang, W. Zhao, T. Zhu, A. Wang, Phys. Rev. D99, 023010 (2019).
  * (48) X. Zhao, C. Zhang, K. Lin, T. Liu, R. Niu, B. Wang, S.-J. Zhang, X. Zhang, W. Zhao, T. Zhu, A. Wang, Gravitational waveforms and radiation powers of the triple system PSR J0337+1715 in modified theories of gravity, Phys. Rev. D100, 083012 (2019).
  * (49) C. Zhang, X. Zhao, A. Wang, B. Wang, K. Yagi, N. Yunes, W. Zhao and T. Zhu, Gravitational waves from the quasicircular inspiral of compact binaries in Einstein-aether theory, Phys. Rev. D101, 044002 (2020).
  * (50) J. Oost, S. Mukohyama, A. Wang, Spherically symmetric exact vacuum solutions in Einstein-aether theory, Universe 7 (2021) 272 [arXiv:2106.09044].
  * (51) T. Gupta, M. Herrero-Valea, D. Blas, E. Barausse, N. Cornish, K. Yagi, N. Yunes, New binary pulsar constraints on Einstein-æther theory after GW170817, Class. Quantum Grav. 38, 195003 (2021).
  * (52) C. Eling and T. Jacobson, Spherical Solutions in Einstein-Aether Theory: Static Aether and Stars, Class. Quantum. Grav. 23, 5625 (2006).
  * (53) C. Eling and T. Jacobson, Black holes in Einstein-aether theory, Class. Quantum. Grav. 23, 5643 (2006).
  * (54) T. Tamaki and U. Miyamoto, Generic features of Einstein-Aether black holes, Phys. Rev. D77, 024026 (2008).
  * (55) D. Blas and S. Sibiryakov, Hořava gravity versus thermodynamics: The black hole case, Phys. Rev. D84, 124043 (2011).
  * (56) E. Barausse, T. Jacobson and T.P. Sotiriou, Black holes in Einstein-aether and Hořava-Lifshitz gravity, Phys. Rev. D83, 124043 (2011).
  * (57) P. Berglund, J. Bhattacharyya and D. Mattingly, Mechanics of universal horizons, Phys. Rev. D85, 124019 (2012).
  * (58) A. Wang, Stationary and slowly rotating spacetimes in Hořava-Lifshitz gravity, Phys. Rev. Lett. 110, 091101 (2013).
  * (59) C. Ding, A. Wang and X. Wang, Charged Einstein-aether black holes and Smarr formula, Phys. Rev. D92, 084055 (2015).
  * (60) C. Zhang, X. Zhao, K. Lin, S.-J. Zhang, W. Zhao and A.-Z. Wang, Spherically symmetric static black holes in Einstein-aether theory, Phys. Rev. D102, 064043 (2020).
  * (61) A. Adam, P. Figueras, T. Jacobson, and T. Wiseman, Rotating black holes in Einstein-aether theory, arXiv:2108.00005.
  * (62) R. Chan, M. F. A. da Silva, V. H. Satheeshkumar, Thermodynamics of Einstein-Aether Black Holes, arXiv:2112.14978.
  * (63) M. Bhattacharjee, S. Mukohyama, M.-B. Wan, and A. Wang, Gravitational collapse and formation of universal horizons in Einstein-$\ae$ther theory, Phys. Rev. D98, 064010 (2018).
  * (64) D. Garfinkle, C. Eling, and T. Jacobson, Numerical simulations of gravitational collapse in Einstein-aether theory, Phys. Rev. D76, 024003 (2007).
  * (65) T. Clifton, P.G. Ferreira, A. Padilla, C. Skordis, Modified Gravity and Cosmology, Phys. Reports 513, 1 (2012).
  * (66) D. Langlois, Dark energy and modified gravity in degenerate higher-order scalar-tensor (DHOST) theories: A review, Inter. J. Mod. Phys. D28, 1942006 (2019).
  * (67) A. Bourgoin, et al., Constraining velocity-dependent Lorentz and CPT violations using lunar laser ranging, Phys. Rev. D103, 064055 (2021).
  * (68) H. Pihan-le Bars, etal., New Test of Lorentz Invariance Using the MICROSCOPE Space Mission, Phys. Rev. Lett. 123, 231102 (2019).
  * (69) C.G. Shao, et al., Combined Search for a Lorentz-Violating Force in Short-Range Gravity Varying as the Inverse Sixth Power of Distance, Phys. Rev. Lett. 122, 011102 (2019).
  * (70) A. Bourgoin, C. Le Poncin-Lafitte, A. Hees, S. Bouquillon, G. Francou, and M.-C. Angonin, Lorentz Symmetry Violations from Matter-Gravity Couplings with Lunar Laser Ranging, Phys. Rev. Lett. 119, 201102 (2017).
  * (71) N.A. Flowers, C. Goodge, and J.D. Tasson, Superconducting-Gravimeter Tests of Local Lorentz Invariance, Phys. Rev. Lett. 119, 201101 (2017).
  * (72) A. Kostelecky and N. Russell, Data tables for Lorentz and CPT violation, Rev. Mod. Phys. 83 11 (2011) [arXiv:0801.0287v15, January 2022 Edition].
  * (73) J. Collins, A. Perez, D. Sudarsky, L. Urrutia, and H. Vucetich, Lorentz invariance and quantum gravity: an additional fine-tuning problem?, Phys. Rev. Lett. 93 (2004) 191301.
  * (74) D. Mattingly, Modern Tests of Lorentz Invariance, Living Rev. Relativity, 8, 5 (2005).
  * (75) S. Liberati, Tests of Lorentz invariance: a 2013 update, Class. Quantum Grav. 30, 133001 (2013).
  * (76) M. Pospelov and C. Tamarit, Lifshitz-sector mediated SUSY breaking, J. High Energy Phys. 01 (2014) 048.
  * (77) A. Wang, Hořava gravity at a Lifshitz point: A progress report, Inter. J. Mod. Phys. D26, 1730014 (2017).
  * (78) J. Bhattacharyya, A. Coates, M. Colombo, and T. P. Sotiriou1, Evolution and spherical collapse in Einstein-Æther theory and Hor̆ava gravity, Phys. Rev. D93, 064056 (2016).
  * (79) J. Bhattacharyya, M. Colombo and T. P. Sotiriou, Causality and black holes in spacetimes with a preferred foliation, Class. Quantum Grav. 33, 235003 (2016).
  * (80) T. Jacobson and D. Mattingly, Einstein-aether waves, Phys. Rev. D70, 024003 (2004).
  * (81) J. W. Elliott, G. D. Moore and H. Stoica, Constraining the New Aether: Gravitational Cherenkov Radiation, JHEP 0508, 066 (2005).
  * (82) J. Calderón Bustillo, P.D. Lasky, and E. Thrane, Black-hole spectroscopy, the no-hair theorem, and GW150914: Kerr versus Occam, Phys. Rev. D103, 024041 (2021).
  * (83) R. Abbott, etal. [The LIGO and Virgo Collaborations], Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog, Phys. Rev. D103, 122002 (2021).
  * (84) X. J. Forteza, S. Bhagwat, P. Pani, and V. Ferrari, Spectroscopy of binary black hole ringdown using overtones and angular modes, Phys. Rev. D102, 044053 (2020).
  * (85) S. Bhagwat, X. J. Forteza, P. Pani, and V. Ferrari, Ringdown overtones, black hole spectroscopy, and no-hair theorem tests, Phys. Rev. D101, 044033 (2020).
  * (86) M. Giesler, M. Isi, M. Scheel, and S. Teukolsky, Black Hole Ringdown: The Importance of Overtones, Phys. Rev. X 9, 041060 (2019).
  * (87) M. Isi, M. Giesler, W. M. Farr, M. A. Scheel, and S. A. Teukolsky,Testing the No-Hair Theorem with GW150914, Phys. Rev. Lett. 123, 111102 (2019).
  * (88) G. Carullo, W. Del Pozzo, and J. Veitch, Observational black hole spectroscopy: A time-domain multimode analysis of GW150914, Phys. Rev. D 99, 123029 (2019).
  * (89) C.D. Capano and A.H. Nitz, Binary black hole spectroscopy: A no-hair test of GW190814 and GW190412, Phys. Rev. D102, 124070 (2020).
  * (90) I. Ota and C. Chirenti, Overtones or higher harmonics? Prospects for testing the no-hair theorem with gravitational wave detections, Phys. Rev. D101, 104005 (2020).
  * (91) F.H. Shaik, J. Lange, S.E. Field, R. O’Shaughnessy, V. Varma, L.E. Kidder, H.P. Pfeiffer, and D. Wysocki, Impact of subdominant modes on the interpretation of gravitational-wave signals from heavy binary black hole systems, Phys. Rev. D101, 124054 (2020).
  * (92) N. Uchikata, T. Narikawa, K. Sakai, H. Takahashi, and Hiroyuki Nakano, Black hole spectroscopy for KAGRA future prospect in O5, Phys. Rev. D102, 024007 (2020).
  * (93) S. Bhagwat, M. Cabero, C.D. Capano, B. Krishnan, and D.A. Brown, Detectability of the subdominant mode in a binary black hole ringdown, Phys. Rev. D102, 024023 (2020).
  * (94) M. Cabero, J. Westerweck, C.D. Capano, S. Kumar, A.B. Nielsen, and B. Krishnan, Black hole spectroscopy in the next decade, Phys. Rev. D101, 064044 (2020).
  * (95) A. Maselli, P. Pani, L. Gualtieri, and E. Berti, Parametrized ringdown spin expansion coefficients: A data-analysis framework for black-hole spectroscopy with multiple events, Phys. Rev. D101, 024043 (2020).
  * (96) T. Islam, A.K. Mehta, A. Ghosh, V. Varma, P. Ajith, and B.S. Sathyaprakash, Testing the no-hair nature of binary black holes using the consistency of multipolar gravitational radiation, Phys. Rev. D101, 024032 (2020).
  * (97) B. P. Abbott, et al., [LIGO/Virgo/KAGRA Collaborations], Tests of General Relativity with GWTC-3, arXiv:2112.06861v1 [gr-qc].
  * (98) S. Tsujikawa, C. Zhang, X. Zhao, A. Wang, Odd-parity stability of black holes in Einstein-Aether gravity, Phys. Rev. D104, 064024 (2021).
  * (99) B. P. Abbott et al. [LIGO Scientific and Virgo Col- laborations], GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017).
  * (100) A. Goldstein et al., An Ordinary Short Gamma-Ray Burst with Extraordinary Implications: Fermi-GBM Detection of GRB 170817A, Astrophys. J. 848, L14 (2017).
  * (101) S. M. Carroll and E. A. Lim, Lorentz-violating vector fields slow the universe down, Phys. Rev. D70, 123525 (2004).
  * (102) J. Müller, J.G. Williams and S.G. Turyshev, Lunar laser ranging contributions to relativity and geodesy, Astrophys. Space Sci. Libr. 349 (2008) 457.
  * (103) J. Greenwald, J. Lenells, J.-X. Lu, V.H. Satheeshkumar, A. Wang, Black holes and global structures of spherical spacetimes in Hořava-Lifshitz theory, Phys. Rev. D84, 084040 (2011). 5
  * (104) D. Blas, O. Pujolas, and S. Sibiryakov, A healthy extension of Hořava gravity, Phys. Rev. Lett. 104, 181302 (2010).
  * (105) D. Blas, O. Pujolas, and S. Sibiryakov, Models of non-relativistic quantum gravity: the good, the bad and the healthy, JHEP 04 (2011) 018.
  * (106) P. Hořava, Membranes at quantum criticality, JHEP 0903 (2009) 020; Quantum gravity at a Lifshitz point, Phys. Rev. D79, 084008 (2009).
  * (107) K. Lin, E. Abdalla, R.-G. Cai and A. Wang, Universal horizons and black holes in gravitational theories with broken Lorentz symmetry, Int. J. Mod. Phys. D23, 1443004 (2014).
  * (108) T. Jacobson, Undoing the twist: The Horava limit of Einstein-aether theory, Phys. Rev. D89, 081501(R) (2014); Phys. Rev. D 81, 101502 (2010); Extended Horava gravity and Einstein-aether theory, Phys. Rev. D81, 101502(R) (2010); D82, 129901(E) (2010).
  * (109) R. M. Wald, General Relativity (University of Chicago Press, 2010).
  * (110) K. Lin, V.H. Satheeshkumar, and A. Wang, Static and rotating universal horizons and black holes in gravitational theories with broken Lorentz invariance, Phys. Rev. D93, 124025 (2016).
  * (111) J. Bhattacharyya, Aspects of holography in Lorentz-violating gravity, PhD Dissertation, https://scholars.unh.edu/cgi/viewcontent.cgi?article
=1728$\&$context=dissertation; E. Barausse1 and T. P. Sotiriou, Black holes in
Lorentz-violating gravity theories, Class. Quantum. Grav. 30, 244010 (2013).

  * (112) K. Lin, O. Goldoni, M.F. da Silva, and A. Wang, A New Look at Those Old Black Holes: Existence of Universal Horizons, Phys. Rev. D91, 024047 (2015).
  * (113) J. E. Thompson, H. Chen and B. F. Whiting, Gauge invariant perturbations of the Schwarzschild spacetime, Class. Quantum Grav. 34 174001 (2017).
  * (114) T. Regge and J. A. Wheeler, Stability of a Schwarzschild Singularity, Phys. Rev. 108, 4 (1957).
  * (115) R. d’ Inverno, Introducing Einstein’s Relativity (Oxford University Press, New York, 2000).
  * (116) E. Butkov, Mathematical Physics (Addison-Wesley Publishing Company, Canada, 1972).
  * (117) C. Zhang, T. Zhu and A. Wang, Gravitational axial perturbations of Schwarzschild-like black holes in dark matter halos, Phys. Rev. D104, 124082 (2021); C. Zhang, T. Zhu, X. Fang and A. Wang, Imprints of dark matter on gravitational ringing of supermassive black holes, Phys. Dark Univ. 37, 101078 (2022).
  * (118) A. Ghosh, R. Brito and A. Buonanno, Constraints on quasinormal-mode frequencies with LIGO-Virgo binary–black-hole observations, Phys. Rev. D103, 124041 (2021).