and imaginary parts of suitably transformed EOM. These components are
invariant under complex conjugation, thus, by taking the real and imaginary
part of the EOM, we are effectively projecting them into a basis of complex
conjugation eigenstates. To generalize the result to Kerr BHs, we followed a
similar approach, projecting the EOM into a basis of eigenstates of the parity
complex conjugation operator $\mathcal{P}$. In both cases, we expected the
approach to work because of the parity invariance of the equations of
motion.141414In the Schwarzschild limit, due to spherical symmetry, when we
take the real and imaginary part of the EOM, we are also indirectly projecting
into a basis of eigenstates of $\mathcal{P}$. We did not make this more
explicit because the EOM look slightly more complicated. In this section, we
are no longer dealing with parity invariant EFT corrections, thus the approach
should break down. Notwithstanding, it may be possible to choose a different
basis that shares the same properties as the ones in section IV.2, where the
equations decouple.

Define $\mathcal{A}_{\delta}$, with $\delta\in[-\pi,\pi)$, as an operator that
takes $\alpha\in\mathbb{C}$ and projects perpendicularly onto the line in the
complex plane defined by $re^{i\delta},r\in\mathbb{R}$. We get:

$\mathcal{A}_{\delta}\,\alpha:=\frac{1}{2}\left(e^{-i\delta}\alpha+e^{i\delta}\bar{\alpha}\right)\,.$
(205)

As expected, $\mathcal{A}_{0}\,\alpha=\mathcal{R}e(\alpha)$ and
$\mathcal{A}_{\frac{\pi}{2}}\,\alpha=\mathcal{I}m(\alpha)$. Furthermore, note
that the output is always a real number, and thus invariant under complex
conjugation. Crucially, if we know the projection of $\alpha$ at two different
angles, $\delta_{1},\delta_{2}$, as long as they do not differ by an integer
multiple of $\pi$, we can always recover $\alpha$ and $\bar{\alpha}$:

$\displaystyle\alpha$
$\displaystyle=2\,\frac{e^{-i\delta_{1}}\mathcal{A}_{\delta_{1}}\,\alpha-e^{-i\delta_{2}}\mathcal{A}_{\delta_{2}}\,\alpha}{e^{-2i\delta_{1}}-e^{-2i\delta_{2}}}\,,$
(206) $\displaystyle\bar{\alpha}$
$\displaystyle=2\,\frac{e^{i\delta_{1}}\mathcal{A}_{\delta_{1}}\,\alpha-e^{i\delta_{2}}\mathcal{A}_{\delta_{2}}\,\alpha}{e^{2i\delta_{1}}-e^{2i\delta_{2}}}\,.$

$\mathcal{A}_{\delta}$, allows for a systematic way to find a basis to express
$(\alpha,\bar{\alpha})$ whose components are real valued. Note that
$\mathcal{A}_{\delta}$ is periodic in $\delta$. For any integer $k$, we have:

$\mathcal{A}_{\delta+k\pi}\,\alpha=e^{ik\pi}\mathcal{A}_{\delta}\,\alpha\,.$
(207)

Furthermore, flipping the sign of $\delta$ is equivalent to taking the complex
conjugate of the argument:

$\mathcal{A}_{-\delta}\,\alpha=\mathcal{A}_{\delta}\,\bar{\alpha}\,.$ (208)

Finally, note that rotating $\alpha$ by some phase $\gamma\in[-\pi,\pi)$, is
the same as rotating the $\delta$ line in the opposite direction, i.e.
$\delta\rightarrow\delta-\gamma$:

$\mathcal{A}_{\delta}e^{i\gamma}\alpha=\mathcal{A}_{\delta-\gamma}\alpha$
(209)

Now, lets employ this approach to a toy problem, similar to equation (131).
Imagine we want to solve the following complex valued differential equation:

$A\psi+\lambda\,B\bar{\psi}=0\,.$ (210)

Here, $A,B$ are real valued differential operators, $\psi$ is a complex valued
function, and $\lambda\in\mathbb{C}$ has phase
$\lambda_{\theta}\in[-\pi,\pi)$:

$\lambda:=\absolutevalue{\lambda}e^{i\lambda_{\theta}}$ (211)

Because the equation depends on both $\psi$ and $\bar{\psi}$, it is not
holomorphic in $\psi$. Thus, we need to explicitly solve the real and
imaginary parts of the equation, and this is in general non trivial. However,
there might be a lines in the complex plane where the equations decouple.
Applying $\mathcal{A}_{\delta}$ to both sides, and using the fact that $A$ and
$B$ commute with complex conjugation, we get:

$A\,\mathcal{A}_{\delta}\,\psi+\absolutevalue{\lambda}B\,\mathcal{A}_{\lambda_{\theta}-\delta}\,\psi=0$
(212)

Now, we want to pick two distinct values of $\delta$ such that (212) depends
on a single real function, i.e.
$\mathcal{A}_{\delta}\psi\sim\mathcal{A}_{\lambda_{\theta}-\delta}$. Using
equation (207), we get two distinct classes of values for $\delta$, labelled
as $\delta_{\pm}$:

$\displaystyle\delta_{+}$
$\displaystyle=-\delta_{+}+\lambda_{\theta}+2k_{+}\pi$ (213)
$\displaystyle\delta_{-}$
$\displaystyle=-\delta_{-}+\lambda_{\theta}+(2k_{-}+1)\pi$

where $k_{\pm}\in\mathbb{Z}$. Solving this, we get:

$\displaystyle\delta_{+}$
$\displaystyle=\frac{\lambda_{\theta}}{2}+k_{+}\pi\,,$ (214)
$\displaystyle\delta_{-}$
$\displaystyle=\frac{\lambda_{\theta}}{2}+\frac{\pi}{2}+k_{-}\pi\,.$

Plugging this into (212) we get two decoupled equations:

$\begin{cases}(A+\absolutevalue{\lambda}B)\,\mathcal{A}_{\delta_{+}}\,\psi=0\\\
(A-\absolutevalue{\lambda}B)\,\mathcal{A}_{\delta_{-}}\,\psi=0\end{cases}$
(215)

We can solve this in terms of $\mathcal{A}_{\delta_{+}}\,\psi$ and
$\mathcal{A}_{\delta_{-}}\,\psi$151515All values of $k_{+},k_{-}$ yield the
same EOM., and then, using equation (206), recover the value of $\psi$ and
$\bar{\psi}$. Note that the second equation can be obtained from the first by
mapping $\absolutevalue{\lambda}\rightarrow-\absolutevalue{\lambda}$.

Working perturbatively in $L/M$, in section IV.3.1, we will use this approach
to correct the QNMs of Schwarzschild BHs. Then, in section IV.3.2, we will
generalize the approach to find the corrections to QNMs of Kerr BHs of any
spin. To do so, we will define an analogue of $\mathcal{A}_{\delta}$ that
depends on $\mathcal{P}$ instead of the complex conjugation operator. In both
cases, we will show that QNM frequencies reduce to the results of section IV.2
with $\lambda_{(e)}$ replaced by
$\absolutevalue{\lambda}=\sqrt{\lambda_{(e)}^{2}+\lambda_{(o)}^{2}}$, in
agreement with equation (120).

####  IV.3.1 Schwarzschild limit

In this section, we will derive the EOM for the corrected QNMs of a
Schwarzschild BHs, when considering mixed parity EFT corrections. In section
IV.2.2, we computed the corrected EOM for $\lambda_{(o)}=0$, non
perturbatively in $L/M$. It is much harder to proceed non-perturbatively for
mixed parity corrections, as there isn’t an obvious way to decouple the EOM.
Nevertheless, validity of the EFT approximation requires $L\ll M$. Thus, only
the leading order correction is physically relevant and we can safely proceed
perturbatively.

Our starting point is equation (134). Making use of the GHP identities that
are specific to spherically symmetric backgrounds (equation (167)), the EOM
simplify. We get:

$\mathcal{M}=\frac{1}{3}\psi_{2}L^{2}\lambda\,J+\mathcal{O}\left(\frac{L}{M}\right)^{4}\,.$
(216)

where:

$\mathcal{M}_{(\alpha)}=\begin{pmatrix}\eth^{\prime}\,\Phi_{0}-\left(\text{\thorn}-2\rho\right)\Phi_{1}\\\
\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)\,\Phi_{1}-\eth\,\Phi_{2}\\\
\left(\text{\thorn}^{\prime}-\rho^{\prime}\right)\,\Phi_{0}-\eth\,\Phi_{1}\\\
\eth^{\prime}\,\Phi_{1}-\left(\text{\thorn}-\rho\right)\,\Phi_{2}\end{pmatrix}_{(\alpha)}\,,$
(217)

and:

$J_{(\alpha)}=\begin{pmatrix}\eth\,\bar{\Phi}_{0}+2\left(\text{\thorn}+\bar{\rho}\right)\bar{\Phi}_{1}\\\
-\eth^{\prime}\,\bar{\Phi}_{2}-2\left(\text{\thorn}^{\prime}+\bar{\rho}^{\prime}\right)\bar{\Phi}_{1}\\\
-\left(\text{\thorn}+2\bar{\rho}\right)\bar{\Phi}_{2}-2\eth\,\bar{\Phi}_{1}\\\
\left(\text{\thorn}^{\prime}+2\bar{\rho}^{\prime}\right)\bar{\Phi}_{0}+2\eth^{\prime}\bar{\Phi}_{1}\end{pmatrix}_{(\alpha)}\,.$
(218)

The EOM are non-holomorphic, as they depend on both $\Phi_{i}$, and
$\bar{\Phi}_{i}$. To tackle this, we will project the equations along suitably
chosen directions of the complex plane, as outlined in the beginning section
IV.3. However, the action of $\mathcal{A}_{\delta}$ breaks GHP covariance.
When we take $\mathcal{A}_{\delta}\Phi_{0}$, we are adding $\Phi_{0}$ and
$\bar{\Phi}_{0}$, GHP scalars with different spin weights. Thus, as was done
in section IV.2.2, we must reduce / increase the spin weight of
$\Phi_{0}/\Phi_{2}$ respectively, to work with spinless scalars. We define the
following set of quantities with $0$ spin weight:

$P_{0}:=\eth^{\prime}\Phi_{0}\,,\quad P_{1}:=\Phi_{1}\,,\quad
P_{2}:=\eth\Phi_{2}\,.$ (219)

Similarly, we must work with EOM of vanishing spin. To set the spin weights of
the $m$ and $\bar{m}$ components of (216) to $0$, we act with $\eth^{\prime}$
and $\eth$ respectively. The resulting EOM will contain factors of
$\eth\eth^{\prime}P_{1}$ and $\eth\eth^{\prime}\bar{P}_{1}$. As argued in
section IV.2.2, there is a relationship between the $\eth\eth^{\prime}$
operator and the spherical harmonics $Y_{\ell m}$, see (172). Using the fact
that in the Schwarzschild limit, $\text{\thorn},\text{\thorn}^{\prime}$, are
exclusively $(t,r)$ derivatives and $\eth,\eth^{\prime}$ act exclusively in
$(x,\phi)$, we can generalize this relation. For $\eta$, a GHP scalar with
spin weight $0$, that can be decomposed as a linear combination of spherical
harmonics with angular quantum number $\ell$:

$\eta(t,r,x,\phi):=\sum_{m=-\ell}^{m=\ell}\eta_{m}(t,r)Y_{\ell m}(x,\phi)\,,$
(220)

we have:

$\eth\eth^{\prime}\eta=-\frac{\ell(\ell+1)}{2r^{2}}\eta\,.$ (221)

Complex conjugation preserves the angular quantum number, thus it is safe to
assume that $P_{1}$ and $\bar{P}_{1}$ obey this condition. In terms of $P_{i}$
the 0 spin-weight EOM are:161616Note that $\widetilde{\mathcal{M}}$ and
$\widetilde{J}$ are no longer the tetrad components of some tensor.

$\widetilde{\mathcal{M}}_{i}\,P_{i}=\frac{1}{3}\psi_{2}L^{2}\lambda\,\widetilde{J}_{i}\bar{P}_{i}+\mathcal{O}\left(\frac{L}{M}\right)^{4}\,.$
(222)

where the summation in $i$ is implicit and:

$\displaystyle\widetilde{\mathcal{M}}_{i}P_{i}$
$\displaystyle:=\begin{pmatrix}-\left(\text{\thorn}-2\rho\right)\,P_{1}+P_{0}\\\\[3.00003pt]
\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)\,P_{1}-P_{2}\\\\[3.00003pt]
\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)\,P_{0}+\frac{\ell(\ell+1)}{2r^{2}}\,P_{1}\\\\[3.00003pt]
-\left(\text{\thorn}-2\rho\right)\,P_{2}-\frac{\ell(\ell+1)}{2r^{2}}\,P_{1}\end{pmatrix}\,,$
(223) $\displaystyle\widetilde{J}_{i}\bar{P}_{i}$
$\displaystyle:=\begin{pmatrix}2\left(\text{\thorn}+\rho\right)\,\bar{P}_{1}\,+\,\bar{P}_{0}\\\\[3.00003pt]
-2\left(\text{\thorn}^{\prime}+\rho^{\prime}\right)\,\bar{P}_{1}-\bar{P}_{2}\\\\[3.00003pt]
-\left(\text{\thorn}+\rho\right)\,\bar{P}_{2}+\frac{\ell(\ell+1)}{r^{2}}\,\bar{P}_{1}\\\\[3.00003pt]
\left(\text{\thorn}^{\prime}+\rho^{\prime}\right)\,\bar{P}_{0}-\frac{\ell(\ell+1)}{r^{2}}\,\bar{P}_{1}\end{pmatrix}\,.$

Note that the LHS depends exclusively on $P_{i}$, whereas the RHS is a
function of $\bar{P}_{i}$. Furthermore, for all $i$,
$\widetilde{\mathcal{M}}_{i}$ and $\widetilde{J}_{i}$ are real operators, the
only complex parameter in (222) is $\lambda$. Thus, the EOM have the same form
as (210), and we can decouple them into two polarizations following the same
procedure.171717Note that in (210) we work with a single scalar $\psi$,
whereas here we have a vector $\psi_{i}$. Notwithstanding, because all
$\widetilde{J}_{i}$ commute with complex conjugation, the same approach works.
We find two phases, $\delta_{\pm}$, that decouple the EOM. Taking
$\lambda_{\theta}$ to be the phase of $\lambda$, we get:

$\displaystyle\delta_{+}$ $\displaystyle=\frac{\lambda_{\theta}}{2}\,,$ (224)
$\displaystyle\delta_{-}$
$\displaystyle=\frac{\lambda_{\theta}}{2}+\frac{\pi}{2}\,.$

The decoupled variables are:

$\mathcal{P}_{i}^{\pm}:=\mathcal{A}_{\delta^{\pm}}P_{i}\,$ (225)

And the EOM separate as:

$\left(\widetilde{\mathcal{M}}_{i}\mp\frac{1}{3}\psi_{2}L^{2}\absolutevalue{\lambda}\,\widetilde{J}_{i}\right)\,\mathcal{P}_{i}^{\pm}=\mathcal{O}\left(\frac{L}{M}\right)^{4}\,.$
(226)

Now, to derive radial decoupled EOM, we follow an identical procedure to the
one in section IV.2.2. First, we factorize the angular and radial
dependence:181818$\mathcal{P}_{i}^{\pm}$ are real variables, but we chose a
complex valued Ansatz. Because the ODE operators are real valued, we can fix
this by taking the real value of the solution in the end.

$\displaystyle\mathcal{P}_{0}^{\pm}(t,r,x,\phi)$ $\displaystyle=e^{-i\omega
t}Y_{\ell m}(x,\phi)\,\rho(r)\,R_{0}^{\pm}(r)\,,$ (227)
$\displaystyle\mathcal{P}_{1}^{\pm}(t,r,x,\phi)$ $\displaystyle=e^{-i\omega
t}Y_{\ell m}(x,\phi)\,R_{1}^{\pm}(r)\,,$
$\displaystyle\mathcal{P}_{2}^{\pm}(t,r,x,\phi)$ $\displaystyle=e^{-i\omega
t}Y_{\ell m}(x,\phi)\,\rho^{\prime}(r)\,R_{2}^{\pm}(r)\,.$

The spherical harmonics factor out, as there are no angular derivatives in
(226). For each polarization, we get a system of four first-order ODEs of
three radial variables. This is over-determined, thus, as was done in IV.2.2,
we can find linear combinations of the first two equations that eliminate
$R_{1}^{\pm}(r)$ and $\differential R_{1}^{\pm}(r)/\differential r$ from the
$3^{rd}$ and $4^{th}$ equations. The remaining two equations can be decoupled
by setting:

$\displaystyle R_{0}^{\pm}(r)$
$\displaystyle=\frac{1}{r\,f(r)}\left(\xi_{\pm}\,u^{\pm}(r)\,+\,\xi_{\mp}\,\,v^{\pm}(r)\right)\,,$
(228) $\displaystyle R_{2}^{\pm}(r)$
$\displaystyle=\frac{1}{r\,f(r)}\left(\xi_{\pm}\,u^{\pm}(r)\,-\,\xi_{\mp}\,\,v^{\pm}(r)\right)\,,$

where:

$\displaystyle\xi_{\pm}$
$\displaystyle:=1\pm\frac{1}{3}\absolutevalue{\lambda}\left(\frac{L}{M}\right)^{2}\left(\frac{M}{r}\right)^{3}\,,$
(229) $\displaystyle f(r)$ $\displaystyle:=1-\frac{2M}{r}\,.$

The EOM reduce to:

$\begin{cases}\derivative{r}v^{\pm}+\Xi_{\pm}\,\frac{i\omega}{f}\,u^{\pm}=\mathcal{O}\left(\frac{L}{M}\right)^{4}\\\
\derivative{r}u^{\pm}+\left(\Xi_{\mp}\frac{i\omega}{f}+\frac{1}{i\omega}V_{\pm}\right)v^{\pm}=\mathcal{O}\left(\frac{L}{M}\right)^{4}\end{cases}\,,$
(230)

where:

$\displaystyle V_{\pm}$
$\displaystyle:=\frac{\ell(\ell+1)}{r^{2}}\left(1\pm\frac{4}{3}\absolutevalue{\lambda}\left(\frac{L}{M}\right)^{2}\left(\frac{M}{r}\right)^{3}\right)\,,$
(231) $\displaystyle\Xi_{\pm}$
$\displaystyle:=1\pm\absolutevalue{\lambda}\frac{2}{3}\left(\frac{L}{M}\right)^{2}\left(\frac{M}{r}\right)^{3}\,.$

Expanding equation (179) to leading order in $\widetilde{\lambda}_{(e)}$, and
replacing $\widetilde{\lambda}_{(e)}$ with
$\left(L/M\right)^{2}\absolutevalue{\lambda}$, we recover (230). Thus, to
leading order in $L/M$, we find agreement between the QNMs obtained directly
for parity even EFT corrections, and the more general mixed parity case. Both
cases are captured by equation (120).

####  IV.3.2 Generalization to Kerr

We can now generalize our conclusions to Kerr BHs of any spin. As before, we
will prove that the corrections to QNM frequencies can be obtained from the
parity even results, section IV.2.3, by simply replacing $\lambda_{(e)}$ with
$\absolutevalue{\lambda}=\sqrt{\lambda_{(e)}^{2}+\lambda_{(o)}^{2}}$. As in
the previous section, the starting point will be equation (134). Following the
procedure outlined in section IV.2 until equation (187), but keeping
$\lambda_{(o)}\neq 0$, we get:

$\begin{cases}\mathcal{O}_{1}\Phi_{0}-\left(\frac{L}{M}\right)^{2}\,\lambda\,\mathcal{D}\psi_{H}=\mathcal{O}\left(\frac{L}{M}\right)^{4}\\\\[5.0pt]
\Phi_{0}-\mathcal{T}_{0}\mathbb{P}\,\mathcal{P}\psi_{H}=\mathcal{O}\left(\frac{L}{M}\right)^{2}\end{cases}\,,$
(232)

where we remind the reader that $\Phi_{0}$ is the component of the corrected
electromagnetic tensor we want to determine, $\psi_{H}$ is a Hertz potential
for the background, $\lambda=\lambda_{(e)}+i\lambda_{(o)}$ is a complex
$\mathcal{O}(1)$ constant, $\mathbb{P}$ is the parity operator, and
$\mathcal{P}$ the conjugate-parity-transform operator (see section II.5).
$\mathcal{D}$ and $\mathcal{T}_{0}$, defined in section IV.2.3, are GHP
operators with spin-weight $s=0$. In the parity even section, we have proved
that $\mathcal{P}$ commutes with all operators in (232):

$[\mathcal{P},\mathcal{O}_{s}]=[\mathcal{P},\mathcal{D}]=[\mathcal{P},\mathcal{T}_{0}]=[\mathcal{P},\mathbb{P}]=0\,.$
(233)

Equation (232) is almost in the same form as equation (210). However, in the
place of $A,B$, differential operators that commute with complex conjugation,
we have operators that commute with $\mathcal{P}$. This motivates the
definition of a new projection operator $\mathcal{B}_{\delta}$, analogous to
$\mathcal{A}_{\delta}$ defined in (212). For some complex function $\eta$, we
have:

$\mathcal{B}_{\delta}\,\eta:=\frac{1}{2}\left(e^{-i\delta}\eta+e^{i\delta}\mathcal{P}\eta\right)\,.$
(234)

As argued at the end of section II.5, $\mathcal{P}$ preserves GHP weight, thus
this is also true for $\mathcal{B}_{\delta}$. Furthermore, note that all the
properties listed for $\mathcal{A}_{\delta}$ at the beginning of section IV.3,
generalize to $\mathcal{B}_{\delta}$, by replacing complex conjugation with
$\mathcal{P}$.

Because $\mathcal{P}$ commutes with all the GHP operators in equation (232),
so does $\mathcal{B}_{\delta}$. Acting with it on both equations, we get:

$\begin{cases}\mathcal{O}_{1}\,\mathcal{B}_{\delta}\Phi_{0}-\left(\frac{L}{M}\right)^{2}\,\absolutevalue{\lambda}\,\mathcal{D}\,\mathcal{B}_{\delta-\lambda_{\theta}}\psi_{H}=\mathcal{O}\left(\frac{L}{M}\right)^{4}\\\\[5.0pt]
\mathcal{B}_{\delta}\Phi_{0}-\mathcal{T}_{0}\mathbb{P}\,\mathcal{B}_{-\delta}\psi_{H}=\mathcal{O}\left(\frac{L}{M}\right)^{2}\end{cases}\,,$
(235)

where $\lambda_{\theta}$ is the phase of $\lambda$. In both equations we see
$\mathcal{B}_{\delta}\Phi_{0}$, however, in the top equation it is related to
$\mathcal{B}_{\delta-\lambda_{\theta}}\psi_{H}$, whereas in the bottom
equation we relate it to $\mathcal{B}_{-\delta}\psi_{H}$. To get decoupled
equations we must have
$\mathcal{B}_{-\delta}\psi_{H}\sim\mathcal{B}_{\delta-\lambda_{\theta}}\psi_{H}$.
This condition is entirely analogous to what we explored in the toy model at
the beginning of section IV.3. We find two decoupled equations, analogous to
(215):

$\begin{cases}\mathcal{O}_{1}\,\mathcal{B}_{\delta^{\pm}}\Phi_{0}\mp\left(\frac{L}{M}\right)^{2}\,\absolutevalue{\lambda}\,\mathcal{D}\,\mathcal{B}_{-\delta^{\pm}}\psi_{H}=\mathcal{O}\left(\frac{L}{M}\right)^{4}\\\\[6.99997pt]
\mathcal{B}_{\delta^{\pm}}\Phi_{0}-\mathcal{T}_{0}\mathbb{P}\,\mathcal{B}_{-\delta^{\pm}}\psi_{H}=\mathcal{O}\left(\frac{L}{M}\right)^{2}\end{cases}$
(236)

with:

$\displaystyle\delta_{+}$
$\displaystyle=\frac{\lambda_{\theta}}{2}+k_{+}\pi\,,$ (237)
$\displaystyle\delta_{-}$
$\displaystyle=\frac{\lambda_{\theta}}{2}+\frac{\pi}{2}+k_{-}\pi\,.$

This reduces to the EOM for the parity even case, equation (192), under the
transformations:191919The factor of $i^{\frac{1\mp 1}{2}}$ is included because
$\Phi_{0}^{\pm}$ and $\psi_{H}^{\pm}$ are parity even / odd scalars, whereas
$\mathcal{B}_{\delta^{\pm}}\Phi_{0}$ and $\mathcal{B}_{-\delta^{\pm}}\psi_{H}$
are both parity even. This is purely technical and has no bearing on the final
result.

$\displaystyle\mathcal{B}_{\delta^{\pm}}\Phi_{0}$ $\displaystyle\rightarrow
i^{\frac{1\mp 1}{2}}\Phi_{0}^{\pm}\,,$ (238)
$\displaystyle\mathcal{B}_{-\delta^{\pm}}\psi_{H}$
$\displaystyle\rightarrow-i^{\frac{1\mp 1}{2}}\psi_{H}^{\pm}\,,$
$\displaystyle\absolutevalue{\lambda}$
$\displaystyle\rightarrow\lambda_{(e)}\,.$

Thus, in agreement with the findings of section IV.3.1, we prove that the
corrections to Kerr QNM frequencies under mixed parity EFT corrections, can be
obtained from the parity even case by taking
$\lambda_{(e)}\rightarrow\sqrt{\lambda_{(e)}^{2}+\lambda_{(o)}^{2}}$. We have
proved that for any Kerr BH background, QNM frequencies are described by
equation (120).

###  IV.4 Results

In this section, we will solve the EOM numerically, to find the corrections to
the QNM frequencies. In section (IV.4.1) we will study the Schwarzschild limit
of parity even EFT corrections, and then, in IV.4.2 we generalize the result
to Kerr BHs of arbitrary spin. As argued in IV.3, at leading order in $L/M$,
the results can be readily generalized to mixed parity corrections by
replacing $\lambda_{(e)}$ with $\sqrt{\lambda_{(e)}^{2}+\lambda_{(o)}^{2}}$

####  IV.4.1 Schwarzschild limit

Figure 9: In this plot, setting $\lambda_{(o)}=0$, we show the corrected
electromagnetic QNM frequency, for Schwarzschild BHs as a function of
$\lambda_{(e)}$. The plot lines show the solutions to equation (179) obtained
non-perturbatively in $L/M$. The plot markers were extracted from figure 4 of
[51]. In that paper, the authors computed the QNMs using a third-order WKB
approximation, valid in the limit $\ell\gg 1$. $\ell=1$ is outside this
regime, thus, as expected, there is very poor agreement between our results
for $\ell=1$. The situation improves for $\ell=2$ and $\ell=3$, but for large
values of $\widetilde{\lambda}_{(e)}$ the difference is still noticeable. As
argued before, the EFT approximation is only valid in the limit $L\ll M$, thus
most of this plot is non-physical.

Using the direct pseudospectral methods outlined in the appendix A, we solved
equation (179) for $\ell=1,2$, $m\in\\{-\ell,\cdots,\ell\\}$, by working non-
perturbatively in $L/M$. We took $\sim 50$ values of
$\widetilde{\lambda}_{(e)}\in[-1,1]$ and obtained the corresponding QNM
frequency for both polarizations, $\omega^{\pm}$. As seen in section IV.2.1,
in the Eikonal limit, for Schwarzschild BHs, the imaginary part of
$\omega_{1}$ vanishes. Thus, for $\ell\neq\infty$, we expect
$\mathcal{I}m(\omega_{1}^{\pm})$ to be small. This implies that, for
Schwarzschild BHs, $\mathcal{I}m(\omega^{\pm})$ will be particularly close to
$\mathcal{I}m(\omega_{0})$. Thus, in this regime, to appropriately estimate
$\omega_{1}^{\pm}$, we must have very precise values of $\omega^{\pm}$ and
$\omega_{0}$. In [51], the authors solved equation (179) for many values of
$\widetilde{\lambda}_{(e)}$, using a third order WKB approximation. This
approach is only valid in the large $\ell$ limit. Although several results in
the literature argue that WKB methods produce reasonable results even for
$\ell=\mathcal{O}(1)$, the imprecision will necessarily grow the closer we get
to $\ell=1$. In figure 9, we compare our results with theirs. Although we find
rough agreement for $\ell=2,3$, for $\ell=1$ we see a large difference. In
fact, the difference between the two results at $\widetilde{\lambda}_{(e)}=0$
is of the same order of magnitude as the maximum correction for
$\widetilde{\lambda}_{(e)}\in[-1,1]$. This is not surprising given the
discussion above.

Parameters | $100M\omega_{1}^{\pm}$
---|---
$\ell$ | Pol. | Part. | Direct | Perturbative | [51]
1 | - | $\mathcal{R}e$ | $1.0567713780$ | $1.0567713780$ | -
$\mathcal{I}m$ | $0.0612235305$ | $0.0612235305$ | $-0.05$
+ | $\mathcal{R}e$ | $-1.0567713780$ | $-1.0567713780$ | -
$\mathcal{I}m$ | $-0.0612235305$ | $-0.0612235305$ | $-0.04$
2 | - | $\mathcal{R}e$ | $1.7486904518$ | $1.7486904518$ | -
$\mathcal{I}m$ | $0.0311341010$ | $0.0311341010$ | $0.03$
+ | $\mathcal{R}e$ | $-1.7486904518$ | $-1.7486904518$ | -
$\mathcal{I}m$ | $-0.0311341010$ | $-0.0311341010$ | $-0.03$
3 | - | $\mathcal{R}e$ | $2.4648927703$ | $-2.4648927703$ | -
$\mathcal{I}m$ | $0.0174061992$ | $0.0174061992$ | $0.02$
+ | $\mathcal{R}e$ | $-2.4648927703$ | $-2.4648927703$ | -
$\mathcal{I}m$ | $-0.0174061992$ | $-0.0174061992$ | $-0.02$
Table 3: In this table, we compare three values for the Schwarzschild limit of
$\omega_{1}$. The Direct column shows the correction estimated using equations
(179) and (239), taking $\widetilde{\lambda}_{(e)}=10^{-9}$. The Perturbative
column contains $\omega_{1}$, computed using the perturbative approach
developed in section IV.2.3, evaluated at $a=0$. The last column has an
estimation of this frequency deduced from the data in figure 4 of [51]. The
results in the first two columns were truncated to $10$ significant digits,
but we actually observed an absolute difference of at most $\sim 10^{-19}$.
Furthermore, as expected, the corrections for the two polarizations are equal
and opposite. As expected, due to the poor accuracy of the WKB approximation
used in [51] at low $\ell$, there is no agreement in the last column for
$\ell=1$. The situation improves with increasing $\ell$.

We have proved, that at leading order in $\widetilde{\lambda}_{(e)}$, the EFT
corrections to the QNM frequencies corresponding to each polarization should
be equal and opposite. The results in figure 9 seem to be consistent with this
notion. To accurately evaluate this, we must specialize to the small
$\widetilde{\lambda}_{(e)}$ regime. Assuming
$\omega^{\pm}\left(\widetilde{\lambda}_{(e)}\right)$ is regular near
$\lambda_{(e)}=0$, we have:

$\displaystyle\omega_{0}$
$\displaystyle=\frac{1}{2}\left(\omega^{\pm}\left(\widetilde{\lambda}_{(e)}\right)+\omega^{\pm}\left(-\widetilde{\lambda}_{(e)}\right)\right)+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\,,$
(239) $\displaystyle\omega_{1}^{\pm}$
$\displaystyle=\frac{1}{2\widetilde{\lambda}_{(e)}}\left(\omega^{\pm}\left(\widetilde{\lambda}_{(e)}\right)-\omega^{\pm}\left(-\widetilde{\lambda}_{(e)}\right)\right)+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\,.$

Plugging in $\widetilde{\lambda}_{(e)}=10^{-9}$ we obtained very precise
values for $\omega_{0}$ and $\omega_{1}^{\pm}$. The background component
describes an electromagnetic QNM of a Schwarzschild BH. This has recently been
calculated with very high precision, using pseudospectral methods in [37]. We
find exact agreement between our results and the frequencies in table 1 of
that paper. In table 3, we show the results for $\omega_{1}^{\pm}$. We find
that $\absolutevalue{\omega_{1}^{+}+\omega_{1}^{-}}\lesssim 10^{-19}$. In the
next section, working perturbatively in $L/M$, we will generalize the results
to Kerr BHs. Setting $a=0$, we get an independent result for
$\omega_{1}^{\pm}$. Once again, we find agreement with at least 19 digits of
precision. We can also estimate the value of $\omega_{1}^{\pm}$ using the data
in [51]. We fit a second-degree polynomial to the data with small
$\widetilde{\lambda}_{(e)}$, and take the derivative at
$\widetilde{\lambda}_{(e)}=0$. Because there isn’t a lot of data in this
region, we get a very rough estimate. This is included in the last column of
3. As in the non-perturbative regime (see figure 9), the results are roughly
consistent with ours for $\ell=2,3$, but there is no agreement at $\ell=1$.

####  IV.4.2 Kerr results

Parameters | $a/M$
---|---
$\ell$ | $m$ | Component | $0$ | $0.25$ | $0.5$ | $0.75$ | $0.95$ | $0.99$ | $0.999$ | $0.9999$ | $0.999975$
$1$ | $-1$ | $M\mathcal{R}e(\omega_{0})$ | $0.248263$ | $0.233724$ | $0.222121$ | $0.212538$ | $0.205897$ | $0.204656$ | $0.204380$ | $0.204352$ | $0.204350$
$M\mathcal{I}m(\omega_{0})$ | $-0.092488$ | $-0.092998$ | $-0.092881$ | $-0.092290$ | $-0.091557$ | $-0.091390$ | $-0.091352$ | $-0.091348$ | $-0.091348$
$100\,M\mathcal{R}e(\omega_{1})$ | $1.056771$ | $0.811775$ | $0.600993$ | $0.451404$ | $0.380628$ | $0.370092$ | $0.367836$ | $0.367612$ | $0.367594$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.061224$ | $0.131237$ | $0.129569$ | $0.079621$ | $0.035260$ | $0.027951$ | $0.026386$ | $0.026231$ | $0.026218$
$1$ | $0$ | $M\mathcal{R}e(\omega_{0})$ | $0.248263$ | $0.249619$ | $0.253933$ | $0.262079$ | $0.272153$ | $0.274311$ | $0.274777$ | $0.274823$ | $0.274827$
$M\mathcal{I}m(\omega_{0})$ | $-0.092488$ | $-0.092009$ | $-0.090329$ | $-0.086330$ | $-0.078557$ | $-0.075949$ | $-0.075305$ | $-0.075240$ | $-0.075234$
$100\,M\mathcal{R}e(\omega_{1})$ | $1.056771$ | $1.051645$ | $1.033499$ | $0.989413$ | $0.915530$ | $0.910598$ | $0.912477$ | $0.912739$ | $0.912762$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.061224$ | $0.044916$ | $-0.009798$ | $-0.131062$ | $-0.376146$ | $-0.473005$ | $-0.497676$ | $-0.500190$ | $-0.500400$
$1$ | $1$ | $M\mathcal{R}e(\omega_{0})$ | $0.248263$ | $0.267305$ | $0.294091$ | $0.337746$ | $0.417719$ | $0.463399$ | $0.489671$ | $0.497136$ | $0.498655$
$M\mathcal{I}m(\omega_{0})$ | $-0.092488$ | $-0.091007$ | $-0.087677$ | $-0.079456$ | $-0.054429$ | $-0.031292$ | $-0.011609$ | $-0.003805$ | $-0.001886$
$100\,M\mathcal{R}e(\omega_{1})$ | $1.056771$ | $1.304920$ | $1.510306$ | $1.547601$ | $0.867720$ | $0.143470$ | $-0.167164$ | $-0.116308$ | $-0.074255$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.061224$ | $-0.087996$ | $-0.307750$ | $-0.521719$ | $-0.236056$ | $0.115480$ | $0.138910$ | $0.061506$ | $0.033062$
$2$ | $-2$ | $M\mathcal{R}e(\omega_{0})$ | $0.457596$ | $0.425879$ | $0.400193$ | $0.378744$ | $0.363875$ | $0.361104$ | $0.360490$ | $0.360428$ | $0.360423$
$M\mathcal{I}m(\omega_{0})$ | $-0.095004$ | $-0.094954$ | $-0.094378$ | $-0.093493$ | $-0.092647$ | $-0.092467$ | $-0.092426$ | $-0.092422$ | $-0.092422$
$100\,M\mathcal{R}e(\omega_{1})$ | $1.748690$ | $1.428612$ | $1.171122$ | $0.970383$ | $0.844319$ | $0.822191$ | $0.817339$ | $0.816857$ | $0.816817$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.031134$ | $0.082368$ | $0.096079$ | $0.087449$ | $0.073683$ | $0.070710$ | $0.070038$ | $0.069971$ | $0.069966$
$2$ | $2$ | $M\mathcal{R}e(\omega_{0})$ | $0.457596$ | $0.498542$ | $0.555399$ | $0.646915$ | $0.813050$ | $0.909991$ | $0.970549$ | $0.990593$ | $0.995284$
$M\mathcal{I}m(\omega_{0})$ | $-0.095004$ | $-0.094089$ | $-0.091131$ | $-0.082512$ | $-0.054668$ | $-0.030053$ | $-0.010655$ | $-0.003482$ | $-0.005264$
$100\,M\mathcal{R}e(\omega_{1})$ | $1.748690$ | $2.132729$ | $2.565175$ | $2.919763$ | $2.176113$ | $0.930436$ | $0.171548$ | $0.036815$ | $0.0170$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.031134$ | $-0.071956$ | $-0.234223$ | $-0.409970$ | $-0.187044$ | $0.057781$ | $0.030259$ | $0.002701$ | $0.0023$
$3$ | $-3$ | $M\mathcal{R}e(\omega_{0})$ | $0.656899$ | $0.608409$ | $0.569188$ | $0.536482$ | $0.513845$ | $0.509631$ | $0.508696$ | $0.508603$ | $0.508595$
$M\mathcal{I}m(\omega_{0})$ | $-0.095616$ | $-0.095418$ | $-0.094724$ | $-0.093755$ | $-0.092865$ | $-0.092679$ | $-0.092636$ | $-0.092632$ | $-0.092632$
$100\,M\mathcal{R}e(\omega_{1})$ | $2.464893$ | $2.045951$ | $1.713706$ | $1.450022$ | $1.278323$ | $1.247517$ | $1.240734$ | $1.240059$ | $1.240002$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.017406$ | $0.067761$ | $0.090180$ | $0.095180$ | $0.092051$ | $0.090966$ | $0.090706$ | $0.090680$ | $0.090678$
$3$ | $3$ | $M\mathcal{R}e(\omega_{0})$ | $0.656899$ | $0.719584$ | $0.806776$ | $0.947457$ | $1.204076$ | $1.355490$ | $1.451973$ | $1.484571$ | $1.492257$
$M\mathcal{I}m(\omega_{0})$ | $-0.095616$ | $-0.094858$ | $-0.092021$ | $-0.083325$ | $-0.054758$ | $-0.029860$ | $-0.010596$ | $-0.003476$ | $-0.001753$
$100\,M\mathcal{R}e(\omega_{1})$ | $2.464893$ | $2.991711$ | $3.640023$ | $4.310268$ | $3.631857$ | $1.918234$ | $0.593560$ | $0.181262$ | $0.089684$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.017406$ | $-0.073354$ | $-0.214206$ | $-0.373040$ | $-0.198956$ | $0.002294$ | $0.009856$ | $0.001438$ | $0.000387$
Table 4: Numerical values for the QNM frequencies and corresponding EFT
correction of electromagnetic modes. All the included digits are significant.

(a)

(b)

(c)

(d)

Figure 10: Here we show the QNM frequency $\ell=m$ of electromagnetic modes,
and corresponding EFT correction for multiple values of $\ell$. Panels a) and
b) contain the real and imaginary parts of the background frequency, whereas
in panels c) and d) we have the corresponding EFT corrections. The solid dark
lines represent the Eikonal limit prediction derived in section IV.2.1. As
expected, with increasing $\ell$, the QNM frequencies approach the large
$\ell$ limit line. Furthermore, in agreement with figures 2 and 12, in the
near extremal limit, $M\omega_{0}\rightarrow\ell/2$ and
$M\omega_{1}\rightarrow 0$.

(a)

(b)

(c)

(d)

Figure 11: In this image we study how quickly the QNM frequencies approach the
Eikonal limit. On the y-axis, we see the logarithmic derivative,
$\gamma=\differential\log\absolutevalue{\omega^{\ell}-\omega^{\infty}}/\differential\log(\ell)$.
Panels $a)$ and $b)$, show the convergence of the real and imaginary part of
the background frequency, whereas $c)$ and $d)$ have the real and imaginary
part of the EFT correction. Because we expect $\mathcal{R}e\omega$ to grow
linearly with $\ell$, we normalized the difference diving by $\ell$. We have
$M\left(\omega^{\ell}-\omega^{\infty}\right)\sim\ell^{-\gamma}$. From equation
(151) we expect that the real part of the background frequency has
$\gamma\rightarrow-1$, while the imaginary part should have
$\gamma\rightarrow-2$. This is consistent with the results in the figure 10.
We did not derive any theoretical results regarding the expected convergence
speed of the EFT correction. Notwithstanding, we find that for both the real
and imaginary parts of $\omega_{1}$, the convergence rate is consistent with
$\ell^{-1}$. Note that there is erratic behaviour of the EFT correction for
low $\ell$ values, particularly for $a/M=3/7$ and $a/M=4/7$. This happens,
because, as seen in figure 10, $\omega_{1}^{\ell}$ crosses
$\omega_{1}^{\infty}$, in this region, taking the logarithm of the difference
to infinity.

(a)

(b)

Figure 12: Here, we show a parametric plot of the lowest lying electromagnetic
QNM frequencies. The left panel shows the background component, while on the
right we have the corresponding EFT correction. The plot markers show the
Schwarzschild limit, obtained by solving the EOM in section IV.2.2. The
parametric lines show the solution to the EOM in section IV.2.3, for BHs with
$a/M\in(0,0.999975)$. We find great agreement between the two results.
Qualitatively, the frequencies share most features of scalar QNMs. In the
Schwarzschild limit, modes with the same $\ell$ share the same frequency and
correction. In the near-extremal limit, for $\ell=m$ modes, the background
frequency approaches $\ell/2$, whereas the correction approaches $0$.

(a)

(b)

Figure 13: Here, we show the near extremal behaviour of $\ell=m$
electromagnetic QNM frequencies. On the left panel, we see the background
components, while the right panel shows the corresponding EFT corrections. As
seen for scalar QNMs, in the limit $a\rightarrow M$, $M\omega_{0}$ approaches
$\ell/2$, while $M\omega_{1}\rightarrow 0$. Contrary to the behaviour at low
spin, in the near-extremal limit, background modes with larger $\ell$ have a
larger imaginary part. In this figure, the lowest lying QNM is the $\ell=m=3$
mode.

To solve equation (203), and obtain the QNM frequencies corresponding
corrections, we employed pseudospectral methods, as described in the appendix
A. The main results are summarized in figure 2 and table 4. We will discuss
qualitatively these results below, but we will first establish their
consistency and accuracy. In figure 4, we computed $\omega_{1}$ for a BH with
$a=0.999975M$ at varying grid sizes. We proved that, with increasing grid
size, the numerical accuracy of $\omega_{1}$ increased exponentially. This is
in line with the expectation for pseudo-spectral methods.

One of our main consistency checks was done in the Schwarzschild limit. In
section IV.2.2, we derived the EOM for the QNMs in this regime, following an
approach that is not perturbative in $\widetilde{\lambda}_{(e)}$. In table 3,
we compared $\omega_{1}$ derived through that approach, with the solutions of
equation (203), with $a=0$. We found perfect agreement within the numerical
accuracy of both approaches. Due to spherical symmetry, in this regime, QNMs
and corresponding corrections are independent of $m$, with the degeneracy
being broken for $a>0$. This can be seen in figures 2 and 12.

On the other hand, our results are consistent with the Eikonal limit
prediction of section IV.2.1. In figures 10 and 14 we show, accordingly, the
$\ell=m$ and $\ell=-m$ QNM frequencies of Kerr BHs for a large range of $\ell$
values, together with the analytic large $\ell$ limit prediction. We find
that, with increasing $\ell$, the numerical curves converge to the analytic
prediction. To quantify the convergence rate, in figures 11 and 15, we used a
logarithmic derivative approach. We find that the real part of the background
frequency converges at a rate of order $\ell^{-1}$, whereas the imaginary part
converges as $\ell^{-2}$. This is consistent with the prediction in [14].
Although we had no results regarding the expected convergence rate for
$\omega_{1}$, we found that the real and imaginary parts approach the Eikonal
limit as $\ell^{-1}$.

Qualitatively, our results are very similar to the ones obtained for the
scalar case (see section III.3). The crucial difference is the splitting of
the QNM frequencies into two families, corresponding to the two possible
polarizations of light. We have:

$\omega^{\pm}=\omega_{0}\pm(-1)^{m+\ell+1}\left(\frac{L}{M}\right)^{2}\sqrt{\lambda_{(e)}^{2}+\lambda_{(o)}^{2}}\,\omega_{1}\,.$
(240)

We see that parity even and parity odd EFTs corrections change the QNM
frequencies in the exact same way. Furthermore, we find that the corrections
corresponding to the two polarizations are equal and opposite. This only
happens because, in the RHS of equation (187), we were able to cancel all
terms proportional to $\Phi_{i}$, leaving only terms involving
$\bar{\Phi}_{i}$. Otherwise, we should expect the corrected frequencies to not
be centered around $\omega_{0}$. This seems to be the case in [33, 34, 35].
Presumably because the EFT corrects the background spacetime, the authors find
that $\omega_{1}^{+}+\omega_{1}^{-}\neq 0$. Notwithstanding, the authors find
that the mixing between parity even and parity odd QNMs is given by the square
root of the sum, in agreement with (240).

Modes with $\ell=m$, also dubbed slowly damped QNMs, have background
frequencies, with the least negative imaginary part. This means they are the
slowest decaying modes, thus dominating the spectrum of the BH ringdown. From
these, for most of the parameter space, the $\ell=m=1$ is the lowest lying,
However, at $a\approx 0.96M$, $\mathcal{I}m(\omega_{0})$ for $\ell=m=2$ and
$\ell=m=3$ cross the value for $\ell=m=1$, and these become the slowest
decaying QNMs.202020From table 4, we see that the $\ell=m=3$ has the largest
imaginary part of the two, thus being the lowest lying QNM. We did not see
this tendency for modes with $\ell>3$, but it could be present sufficiently
close to extremality. Coincidentally, except very close to extremality, for
each value $\ell$, $\omega_{1}$ has the largest absolute real and imaginary
part when $\ell=m$. These hit a maximum for $a\sim 0.75M$. Thus. these modes
seem to be the most relevant if we wish to experimentally find an upper bound
for $\lambda$.

In figure 12, we show the data of figure 2 in a parametric form, as a function
of $a/M$. The plot markers represent QNM frequencies in the Schwarzschild
limit, as obtained in section IV.4.1, whereas the plot lines correspond to the
QNM frequency when $a\neq 0$. There is a perfect agreement between the two.
Particularly relevant, is the shape of the slowly damped QNMs. Their
parametric lines cover a large region of the complex plane, varying
particularly quickly in the near extremal lime. For these modes, should we
exclude $a\gtrsim 0.99M$, we would lose approximately the last tenth of the
parametric lines. Thus, in figure 13, we zoom in this region of parameter
space, computing the QNMs for BHs with $a/M\in(0,99,0.999975)$. As predicted
in [38], as $a\rightarrow M$, the background QNM frequencies seem to converge
polynomially to $\ell/2$. Analogously, in the same limit, we find that
$M\omega_{1}$ approaches $0$ at a similar convergence rate.212121In the
$\ell=m=1$ and $\ell=m=2$ cases, modes seem to converge slightly quicker than
expected, but, for $\ell=m=3$, the convergence rate is consistent with the
expectation: $\omega_{0}^{\ell=m=3}-\ell/2\sim\sqrt{M-a}$. We think that, by
getting even closer to extremality, we should get the correct convergence
rate. It was harder to quantify the convergence rate of $\omega_{1}$ as a
logarithm derivative approach did not seem to converge for $a\lesssim
0.999975M$. Furthermore, we find that the absolute value of
$M\mathcal{R}e\omega_{1}$ grows with $\ell$. This is consistent with the
prediction of [38] for the background modes. Finally, note that at $a\sim
0.992$, $M\mathcal{R}e\,\omega_{1}$ crosses $0$.

##  V Discussion

In this paper, we computed the correction to BH QNMs, under scalar and
electromagnetic, leading order EFT corrections. In the electromagnetic case,
we proved that the leading order parity even and parity odd EFT corrections
have the same effect on QNM frequencies. Taking $\lambda{(e)/(o)}$ to be the
coupling parameter that control parity even / odd corrections, the combined
effect of the two is controlled by
$\sqrt{\lambda_{(e)}^{2}+\lambda_{(o)}^{2}}$. Furthermore, we found that the
degeneracy between the two polarizations of light is broken, splitting the
QNMs into two families. At leading order, the corrections to the corresponding
frequencies were shown to be equal and opposite.

Electromagnetic QNMs from a BH merger have not yet been detected. Thus, from
an observational perspective, our results are not very relevant. However, as
seen in [33, 34, 35], the main features generalize to pure gravity EFTs
corrections. There, the degeneracy breaking between the two possible
polarizations is particularly relevant. In theory, if we know the precise
direction of a GW source, combining the signal from several GW detectors, we
may disentangle the two GW tensor polarizations, see [68]. Thus, if both
families of QNMs are excited in a BH merger, we may be able to disentangle the
QNM spectrum into contributions corresponding to each polarization, and use
the difference to infer an upper bound for the size of the EFT correction.

We used, scalar, and electromagnetic QNMs as toy models to learn how to, in
general, compute the EFT corrections to gravitational QNMs of Kerr BHs. There,
the challenge is that not only the EOM are corrected, but the background
stationary metric is also corrected. [33, 34, 35] have overcome this
challenge, by working perturbatively in the BH spin. In future work, we expect
that it is possible to extend our method to compute the gravitational QNMs of
BHs with arbitrary spin. The correction to the background metric can be
determined numerically and then, the resulting EOM can be decoupled and solved
with an extension of both our methods.

##  VI Acknowledgements

The author thanks the contribution of his supervisor Harvey S. Reall, in
motivating, discussing, and reviewing the contents of this paper, Andy Zhao,
for initial calculations regarding the correction to electromagnetic
quasinormal modes in the Schwarzschild limit, and Jorge E. Santos for helpful
discussions regarding the numerical approach. The author also thanks David
Baker, Rifath Khan, Frank Qu, Gonçalo Regado, and Bilyana Tomova for fruitful
physics discussions. The author is jointly funded by the IOA (University of
Cambridge), the University of Cambridge Trust, and King’s College (University
of Cambridge).

##  Appendix A Numerical methods

(a)

(b)

(c)

(d)

Figure 14: Here we show the QNM frequencies of $\ell=-m$ electromagnetic modes
for several values of $\ell$. In panels $a)$ and $b)$, we have the real and
imaginary part of the background frequency, with panels $c)$ and $d)$ showing
the corresponding correction. The thick dark line shows the Eikonal limit
prediction. As expected, the QNM frequencies converge to the asymptotic
prediction.

In this appendix, we will describe the numerical methods used throughout the
paper. We will be using collocation methods, as described in [17, 18, 19, 21].
Effectively, these map continuous ODEs and PDEs problems into finite-
dimensional linear algebra problems.

The final EOM for scalar QNMs in section III and electromagnetic QNMs in
section IV can be written as:

$\displaystyle L^{(0)}(\omega_{0})\,\Phi^{(0)}=0\,,$ (241a) $\displaystyle
L^{(0)}(\omega_{0})\,\Phi^{(1)}+\omega_{1}\delta
L\Phi^{(0)}=L^{(1)}\Phi^{(0)}\,.$ (241b)

As seen in equation (48), $L^{(0)}$ separates as the sum of two ODE operators
and, as such, $\Phi^{(0)}(r,x)=R(r)\Theta(x)$:

$L^{(0)}\Phi^{(0)}(r,x)=\Theta(x)L^{(0)}_{r}R(r)+R(r)L^{(0)}_{x}\Theta(x)\,.$
(242)

QNM frequencies are values of $\omega_{0/1}$ such that there are solutions
$\Phi^{(0/1)}$ that are smooth in the PDE’s domain. In section A.1, we will
explain the main features of the numerical methods used. Then, in section A.2,
we explain how to use this to solve equation (241a), and then, in section A.3,
we will solve equation (241b).

(a)

(b)

(c)

(d)

Figure 15: In this image, we study how quickly do $\ell=-m$ QNM frequencies
converge to the Eikonal limit prediction. As in figure 11, we find that the
real and imaginary parts of the background converge as $\ell^{-1}$, and
$\ell^{-2}$. This is consistent with the prediction in equation (151).
Furthermore, we find that the real and imaginary of the correction converge to
the theoretical prediction as $\ell^{-1}$.

###  A.1 Pseudo-spectral methods

Take a set of functions $\\{T_{n}\\}_{n=1}^{N}$ that are orthonormal with
respect to some inner product defined in a given domain. A common approach to
solving ODE or PDE equations, is to expand its solution in this functional
basis. Then, by understanding the behaviour of differential operators on the
basis of functions, we can map the continuous differential equation onto a
discrete problem. In general, we expect that by increasing $N$, higher
frequency features will be captured and the discrete solution converges to the
true continuous solution. In the literature, this approach is usually known as
the Galerkin, or spectral method. Another possible approach is to sample the
differential equation’s solution in a discrete collocation grid
$\\{x_{n}\\}_{n=1}^{N-1}$, that spans its domain. Then, we can map the
differential operators into matrix operators, reducing the continuous problem
into a simpler finite-dimensional version. If $\\{T_{n}\\}_{n=1}^{N}$ is a set
of polynomials in ascending order of degree, there is a correspondence between
the two methods. Choosing $\\{x_{n}\\}_{n=1}^{N-1}$ to be the $N-1$ roots of
$\\{T_{n}\\}_{n=1}^{N}$, the solutions derived from the Galerkin and
collocation methods are identical. Hence, collocation methods are many times
known in the literature as pseudo-spectral methods.

It has been shown that, in general, that when $\phi_{n}$ are the Chebyshev
polynomials, pseudo-spectral methods have exponential convergence properties.
Furthermore, because the corresponding grid heavily samples values close to
the boundaries of the domain, it is especially suited to BVPs, see [21].
Throughout this paper, we will work with this basis. $T_{n}$ will denote the
$n^{th}$ Chebyshev polynomial, whereas $\\{x_{i}\\}$ encodes the corresponding
Chebyshev grid:

$x_{i}=\frac{1}{2}\left(1-\cos\left(\frac{i\pi}{N}\right)\right),\qquad
i=0,\cdots,N\,.$ (243)

Here, the $N-1$ values in the middle of the grid are the roots of $T_{N-2}$,
whereas $x_{0}$ and $x_{N}$ are needed to enforce the BCs of our problem.

Lets first understand how to use this approach to discretize ODEs. Take:

$L(x)f(x)=0\,,$ (244)

where

$L(x):=\sum_{n}L^{(n)}(x)\derivative[n]{x}\,,$ (245)

is some differential operator, containing functions of $x$ and total
derivatives with respect to $x$. To obtain the corresponding discrete version,
we start by choosing some Chebyshev grid $x_{i}$. Then, sample $f(x)$ in the
Chebyshev grid, defining $f_{i}:=f(x_{i})$. Differentiation is a linear
operation, thus, it must have a matrix representation in the discrete space.
There must be some operator $D$ such that:

$f^{\prime}(x_{i})=\sum_{j}D_{ij}f_{j}\,.$ (246)

In fact, for a Chebyshev grid, the expression of $D$ can be found in section 6
of [20]. The discretization of $L(x)$ can then be obtained by replacing
$\derivative[n]{x}$ with $D^{n}$ and $L_{n}(x)$ with diagonal matrices
$\hat{L}^{(n)}_{ij}:=L^{(n)}(x_{i})\mathbb{I}_{ij}$, where $\mathbb{I}$ is the
identity matrix. Then, the discrete representation is:

$\sum_{j}\hat{L}_{ij}f_{j}:=\sum_{n,k,j}\hat{L}^{(n)}_{ik}(D^{n})_{kj}f_{j}\,.$
(247)

Treating $f_{i}$ as an unknown, we can solve the equation directly using
standard methods.

The generalization to PDEs only requires small changes. Take

$L(x^{1},\cdots,x^{n})\,\Phi(x^{1},\cdots,x^{n})=0\,,$ (248)

where $L$ is a PDE operator with $n$ variables. Start by defining $n$
Chebyshev grids $\\{x^{k}_{i}\\}$, that need not have the same size. The
discretization of $\Phi$ will be a $n$-tensor

$\Phi_{i_{1},\cdots,i_{n}}:=\Phi(x_{i_{1}}^{1},\cdots,x_{i_{n}}^{n})\,.$ (249)

Then, the discrete version of partial derivative operators will be:

$\partialderivative{x^{k}}\rightarrow\underbrace{\mathbb{I}\otimes\cdots\otimes\mathbb{I}}_{k-1}\otimes
D\otimes\underbrace{\mathbb{I}\otimes\cdots\otimes\mathbb{I}}_{n-k}\,.$ (250)

Mapping the remaining factors to ”diagonal tensors”, the discretization of
(248) will take the form:

$\hat{L}_{i_{1},\cdots,i_{n}}^{j_{1},\cdots,j_{n}}\Phi_{j_{1},\cdots,j_{n}}=0$
(251)

where the Einstein summation convention was assumed.

When using pseudo-spectral methods, it is important to quantify the size of
the numerical error. Because these methods converge exponentially with the
grid size, we can simply solve the EOM using two different grid sizes. The
absolute difference between the values obtained for each is a good estimate of
the error size. In this paper, in most cases, we set the second grid to be
roughly $10\%$ bigger than the first.

###  A.2 Solving the background equation

Following the discretization procedure, we can now solve the EOM and find the
QNMs. We will first describe how to solve the background equation, equation
(241a). As argued before, this equation separates into a system of two ODEs.
Choosing a Chebyshev grid with size $N_{x}$ in the angular direction and a
grid with size $N_{r}$ in the radial direction, the discrete representation of
equation (241a) is:

$\displaystyle A_{x}(J,\,m,\,\mathcal{A}_{\ell
m},\,\omega_{0})\cdot\underline{\Theta}=0\,,$ (252a) $\displaystyle
A_{r}(J,\,m,\,\mathcal{A}_{\ell m},\,\omega_{0})\cdot\underline{R}=0\,.$
(252b)

Here $A_{x}$ and $A_{r}$ are matrix operators, $J$ and $m$ are parameters of
the problem (see equations (14) and (46)), $\mathcal{A}_{\ell m}$ and
$\omega_{0}$ are to be treated as eigenvalues, and $\underline{R}$ and
$\underline{\Theta}$ are the corresponding eigenvectors. Throughout the rest
of this section, all underlined quantities are to be understood as vectors.
The two equations depend on both eigenvalues, hence the system is not
completely decoupled. Added difficulty comes from the fact that the second
equation is quadratic in $\omega$.

To alleviate this, we can start by solving the EOM in the Schwarzschild limit.
As argued in section II.3, in this limit, the angular equation can be solved
analytically, yielding the spin-weighted spherical harmonics.
$\mathcal{A}_{\ell m}$ reduces to $l(l+1)-s(s+1)$. The symmetry of the problem
also removes the $m$ dependence from the second equation. Following section
III.C of [19], we used a linearization approach to reduce (252b) to a standard
eigenvalue problem. The resulting discrete eigenvalue problem can be inserted
in a numeric eigenvalue solver. We used Mathematica’s Eigensystem method,
selecting the Arnoldi method. 222222Most eigenvalues obtained in this manner
are not actual eigenvalues of the continuous system, to find the true
frequencies the computation must be repeated with multiple grid sizes, keeping
only the ones that are consistent throughout. To achieve accurate results, the
use of Mathematica’s arbitrary precision capabilities was crucial. In general,
we used at least 200 digits of precision.

To tackle the case with arbitrary spin, we used the Newton-Raphson approach,
(see section VI.A of [19]). This is an iterative method to find the roots of
any map $f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$. Define
$\underline{x}_{0}$ such that:

$f(\underline{x}_{0})=0$ (253)

Taylor expanding $f(\underline{x})$ around $\underline{x}_{0}$, we get:

$f(\underline{x})=f(\underline{x}_{0})+\delta
f(\underline{x})\cdot\delta\underline{x}+\mathcal{O}\left(\delta\underline{x}^{2}\right),$
(254)

where $\delta\underline{x}:=(\underline{x}-\underline{x}_{0})$ and $\delta
f(\underline{x}):=\partialderivative{f(\underline{x})}{\underline{x}}$. By
definition $f(\underline{x}_{0})=0$, thus, we can approximate
$\delta\underline{x}$ with:

$\delta\underline{x}=\left(\delta
f(\underline{x})\right)^{-1}\,\cdot\,f(\underline{x})\,.$ (255)

We update $\underline{x}\rightarrow\underline{x}-\delta\underline{x}$, and
find a new estimate of $\delta\underline{x}$, iterating the process until
convergence. As long as the initial seed is close enough to
$\underline{x}_{0}$, this method has quadratic convergence. However, if the
seed is too distant from $\underline{x}_{0}$, the method will either not
converge, or converge to a different root of $f$.

To apply this method to (252b), we must first fix some leftover gauge freedom.
The vectors $\underline{\Theta}$ and $\underline{R}$ can be rescaled by
constant factors, thus, we must supplement the EOM with two normalization
conditions. We choose two vectors $\underline{n}^{(x)}$ and
$\underline{n}^{(r)}$, and require:

$\begin{cases}\underline{n}^{(x)}\,\cdot\,\underline{\Theta}=1\,,\\\
\underline{n}^{(r)}\,\cdot\,\underline{R}=1\,.\end{cases}$ (256)

Taking
$\underline{x}=\left(\underline{\Theta},\underline{R},\omega_{0},\mathcal{A}_{\ell
m}\right)$, we use equations (252b) and (256) to define $f$. The Newton-
Raphson update equation will be:

$\begin{pmatrix}A_{x}&\partialderivative{\omega_{0}}A_{x}\,\cdot\,\underline{\Theta}&\partialderivative{\mathcal{A}_{\ell
m}}A_{x}\,\cdot\,\underline{\Theta}\\\
A_{r}&\partialderivative{\omega_{0}}A_{r}\,\cdot\,\underline{R}&\partialderivative{\mathcal{A}_{\ell
m}}A_{r}\,\cdot\,\underline{R}\\\ \underline{n}^{(x)}&0&0\\\
\underline{n}^{(r)}&0&0\end{pmatrix}\cdot\begin{pmatrix}\underline{\delta\Theta}^{T}\\\
\underline{\delta R}^{T}\\\ \delta\omega_{0}\\\ \delta\mathcal{A}_{\ell
m}\end{pmatrix}\,=\,\begin{pmatrix}\left(A_{x}\cdot\underline{\Theta}\right)^{T}\\\
\left(A_{r}\cdot\underline{R}\right)^{T}\\\
\underline{n}^{(x)}\cdot\underline{\Theta}-1\\\
\underline{n}^{(r)}\cdot\underline{R}-1\end{pmatrix}$ (257)

Now, to find the background QNM frequency for a background of any spin, we
just need an accurate enough initial guess. We take the Schwarzschild limit
frequency and use it as a seed for a very slowly rotating BH. After
convergence, we can use the resulting frequency to obtain the QNM of a BH with
a slightly larger spin. We iterate this, increasing $J$ gradually, until we
cover the whole range of Kerr BHs, $J\in[0,1)$.

Deciding how much to increase $J$ at each iterative step, is subtle. To aid in
this decision, and speed-up convergence, we can estimate
$\differential\underline{x}_{0}/\differential J$. Including $J$ in equation
(253), we have $f(\underline{x}_{0},J)=0$. Thus, we can treat
$\underline{x}_{0}$ as a function of $J$, parametrizing the level curves of
$f$. Along these curves, we have:

$0=\derivative{J}f\left(\underline{x}_{0}(J),J\right)=\partialderivative{f}{\underline{x}_{0}}\,\cdot\,\derivative{x_{0}}{J}+\partialderivative{f}{J}\,.$
(258)

Solving with respect to $\differential\underline{x}_{0}/\differential J$ we
get:

$\derivative{\underline{x}_{0}}{J}=-\left(\delta
f\right)^{-1}\,\cdot\,\partialderivative{f}{J}\,.$ (259)

Plugging this in the definition of the derivative, we can estimate
$\underline{x}_{0}(J+\delta J)$:

$\underline{x}_{0}(J+\delta J)\simeq-\underline{x}_{0}(J)-\left(\delta
f\right)^{-1}\,\cdot\,\partialderivative{f}{J}\,\delta J\,.$ (260)

In our calculations, we fixed the size of $\delta\mathcal{I}m(\omega_{0})$,
deducing the corresponding size of $\delta J$.

###  A.3 Solving the perturbation equation

Using the background QNMs, we can now solve the perturbation equation,
equation (241b). Following the prescription in section A.1 we obtain the
corresponding discrete version. Because the EOM is a 2 dimensional PDE,
$\Phi^{(0/1)}$ are represented by $N_{x}\crossproduct N_{r}$ matrices, and
$L^{(0)}$, $L^{(1)}$ and $\delta L$ will be represented by 4-tensors with size
$N_{x}\crossproduct N_{r}\crossproduct N_{x}\crossproduct N_{r}$. Fixing
normalization freedom as was done (256), we get:

$\displaystyle L^{(0)}_{ijkl}\Phi_{kl}^{(1)}+\omega_{1}\delta
L^{(1)}_{ijkl}\Phi_{kl}^{(0)}=\mathcal{L}^{(1)}_{ijkl}\Phi_{kl}^{(0)},5em.$
(261a)
$\displaystyle\underline{n}^{(x)}_{i}\cdot\left(\Phi^{(0)}_{ij}+\lambda\Phi^{(1)}_{ij}\right)\cdot\underline{n}^{(r)}_{j}=1$
(261b)

where we assumed the Einstein summation convention. These tensors depend on
the same parameters as the background and on $(\omega_{0},\mathcal{A}_{\ell
m})$. Furthermore, the separability condition at $0$ order translates to
$\Phi^{(0)}=\underline{\Theta}\otimes\underline{R}$. Using equation (256), the
normalization condition (261b) simplifies to
$\underline{n}^{(x)}\cdot\Phi^{(1)}\cdot\underline{n}^{(r)}=0$. Using
Mathematica’s Flattening operator, we get a vector representation of
$\Phi^{(0/1)}$. This in turn implies a matrix representation of the 4-tensors
above. In this basis, the discrete EOM take the form of a linear system:

$\mathcal{M}\cdot\underline{x}=\underline{b}\,,$ (262)

where:

$\displaystyle\underline{x}=\left(\underline{\Phi^{(1)}},\,\omega_{1}\right)\,,$
(263) $\displaystyle\mathcal{M}=\begin{pmatrix}L^{(0)}&\underline{\left(\delta
L^{(1)}\cdot\Phi^{(0)}\right)}^{T}\\\\[5.0pt] \underline{n^{(r)}\otimes
n^{(x)}}&0\end{pmatrix}\,,$
$\displaystyle\underline{b}^{T}=\begin{pmatrix}\underline{\left(L^{(1)}\cdot\Phi^{(0)}\right)}^{T}\\\\[5.0pt]
0\end{pmatrix}\,.$

Underlined quantities are $N_{x}\cdot N_{r}$ vectors whereas not underlined
quantities are scalars or $(N_{x}\cdot N_{r})\crossproduct(N_{x}\cdot N_{r})$
matrices, depending on the context. The system (262) can then directly be
solved with Mathematica’s LinearSolve method. In general $\mathcal{M}$ is non-
singular, thus, we expect that there will be a single $\omega_{1}$
corresponding to each background frequency $\omega_{0}$.

Chebyshev differentiation matrices are very dense, thus, in general, the
discrete systems we obtain have very high condition number. Thus, to solve the
EOM we must work with very high-precision arithmetic. This is particularly
relevant in the near extremal limit, where the condition number increases very
fast. $\mathcal{M}$ and $b$ are functions of $\omega_{0}$ and
$\mathcal{A}_{\ell m}$, thus the system’s precision is determined by the
precision of these eigenvalues. To accurately solve the EOM in the near
extremal regime, we needed these eigenvalues to have at least 70 digits of
precision. To achieve such feat, we had to solve the background EOM using very
large Chebyshev grids. In the most extreme scenario, we used $N_{x}=250$ and
$N_{r}=650$. When solving the background EOM, this implies that each step of
the Newton Raphson method involves solving a $\sim 900\crossproduct 900$
system of equations, which is feasible with current computing capabilities.
However, when discretizing the first order PDE, rank of $\mathcal{M}$ will be
$N_{r}\cdot N_{x}\sim 160000$. Even though this matrix is not very dense,
solving this system would be very challenging. To work around this issue, we
can use a clever form of interpolation to reduce the rank of $\mathcal{M}$. As
mentioned in section A.1, there is a correspondence a Chebyshev collocation
grid and the spectral decomposition in a Chebyshev polynomial basis. In fact,
these coefficients are related by a linear map that can be found in [21]. To
project the perturbation equations into a smaller Chebyshev grid, we find
their spectral representation, truncate the spectral sum at a lower order, and
then, map it back to a Chebyshev collocation grid.

## References

  * Regge and Wheeler [1957] T. Regge and J. A. Wheeler, Stability of a Schwarzschild singularity, Phys. Rev. 108, 1063 (1957).
  * Zerilli [1970] F. J. Zerilli, Effective potential for even parity Regge-Wheeler gravitational perturbation equations, Phys. Rev. Lett. 24, 737 (1970).
  * Teukolsky [1973] S. A. Teukolsky, Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electromagnetic and neutrino field perturbations, Astrophys. J. 185, 635 (1973).
  * Bardeen and Press [1973] J. M. Bardeen and W. H. Press, Radiation fields in the schwarzschild background, J. Math. Phys. 14, 7 (1973).
  * Chandrasekhar and Detweiler [1975] S. Chandrasekhar and S. L. Detweiler, The quasi-normal modes of the Schwarzschild black hole, Proc. Roy. Soc. Lond. A 344, 441 (1975).
  * Glampedakis  _et al._ [2017] K. Glampedakis, A. D. Johnson, and D. Kennefick, Darboux transformation in black hole perturbation theory, Phys. Rev. D 96, 024036 (2017), arXiv:1702.06459 [gr-qc] .
  * Vishveshwara [1970] C. V. Vishveshwara, Scattering of Gravitational Radiation by a Schwarzschild Black-hole, Nature 227, 936 (1970).
  * Press [1971] W. H. Press, Long Wave Trains of Gravitational Waves from a Vibrating Black Hole, Astrophys. J. Lett. 170, L105 (1971).
  * Nollert [1999] H.-P. Nollert, TOPICAL REVIEW: Quasinormal modes: the characteristic ‘sound’ of black holes and neutron stars, Class. Quant. Grav. 16, R159 (1999).
  * Kokkotas and Schmidt [1999] K. D. Kokkotas and B. G. Schmidt, Quasinormal modes of stars and black holes, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, and A. O. Starinets, Quasinormal modes of black holes and black branes, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Konoplya and Zhidenko [2011] R. A. Konoplya and A. Zhidenko, Quasinormal modes of black holes: From astrophysics to string theory, Rev. Mod. Phys. 83, 793 (2011), arXiv:1102.4014 [gr-qc] .
  * Mashhoon [1985] B. Mashhoon, Stability of charged rotating black holes in the eikonal approximation, Phys. Rev. D 31, 290 (1985).
  * Dolan [2010] S. R. Dolan, The Quasinormal Mode Spectrum of a Kerr Black Hole in the Eikonal Limit, Phys. Rev. D 82, 104003 (2010), arXiv:1007.5097 [gr-qc] .
  * Yang  _et al._ [2012] H. Yang, D. A. Nichols, F. Zhang, A. Zimmerman, Z. Zhang, and Y. Chen, Quasinormal-mode spectrum of Kerr black holes and its geometric interpretation, Phys. Rev. D 86, 104006 (2012), arXiv:1207.4253 [gr-qc] .
  * Leaver [1985] E. W. Leaver, An Analytic representation for the quasi normal modes of Kerr black holes, Proc. Roy. Soc. Lond. A 402, 285 (1985).
  * Dias  _et al._ [2010] O. J. C. Dias, P. Figueras, R. Monteiro, H. S. Reall, and J. E. Santos, An instability of higher-dimensional rotating black holes, JHEP 05, 076, arXiv:1001.4527 [hep-th] .
  * Miguel [2021] F. S. Miguel, Interior Quasinormal Modes and Strong Cosmic Censorship, Phys. Rev. D 103, 064077 (2021), arXiv:2012.10455 [gr-qc] .
  * Dias  _et al._ [2016] O. J. C. Dias, J. E. Santos, and B. Way, Numerical Methods for Finding Stationary Gravitational Solutions, Class. Quant. Grav. 33, 133001 (2016), arXiv:1510.02804 [hep-th] .
  * Trefethen [2000] L. N. Trefethen,  _Spectral Methods in MATLAB_ (Society for Industrial and Applied Mathematics, 2000) https://epubs.siam.org/doi/pdf/10.1137/1.9780898719598 .
  * Boyd [2000] J. P. Boyd,  _Chebyshev and Fourier Spectral Methods_ , 2nd ed. (Dover Publications, Inc., Mineola, New York, USA, 2000) p. 594.
  * Baibhav  _et al._ [2023] V. Baibhav, M. H.-Y. Cheung, E. Berti, V. Cardoso, G. Carullo, R. Cotesta, W. Del Pozzo, and F. Duque, Agnostic black hole spectroscopy: quasinormal mode content of numerical relativity waveforms and limits of validity of linear perturbation theory, (2023), arXiv:2302.03050 [gr-qc] .
  * Dreyer  _et al._ [2004] O. Dreyer, B. J. Kelly, B. Krishnan, L. S. Finn, D. Garrison, and R. Lopez-Aleman, Black hole spectroscopy: Testing general relativity through gravitational wave observations, Class. Quant. Grav. 21, 787 (2004), arXiv:gr-qc/0309007 .
  * Brito  _et al._ [2018] R. Brito, A. Buonanno, and V. Raymond, Black-hole Spectroscopy by Making Full Use of Gravitational-Wave Modeling, Phys. Rev. D 98, 084038 (2018), arXiv:1805.00293 [gr-qc] .
  * Ghosh  _et al._ [2021] A. Ghosh, R. Brito, and A. Buonanno, Constraints on quasinormal-mode frequencies with LIGO-Virgo binary–black-hole observations, Phys. Rev. D 103, 124041 (2021), arXiv:2104.01906 [gr-qc] .
  * Silva  _et al._ [2022] H. O. Silva, A. Ghosh, and A. Buonanno, Black-hole ringdown as a probe of higher-curvature gravity theories, (2022), arXiv:2205.05132 [gr-qc] .
  * Weinberg [2013] S. Weinberg,  _The quantum theory of fields. Vol. 2_ (Cambridge University Press, 2013).
  * Burgess [2004] C. P. Burgess, Quantum gravity in everyday life: General relativity as an effective field theory, Living Rev. Rel. 7, 5 (2004), arXiv:gr-qc/0311082 .
  * Burgess [2007] C. P. Burgess, Introduction to Effective Field Theory, Ann. Rev. Nucl. Part. Sci. 57, 329 (2007), arXiv:hep-th/0701053 .
  * Reall and Warnick [2022] H. S. Reall and C. M. Warnick, Effective field theory and classical equations of motion, J. Math. Phys. 63, 042901 (2022), arXiv:2105.12028 [hep-th] .
  * Cardoso  _et al._ [2018] V. Cardoso, M. Kimura, A. Maselli, and L. Senatore, Black Holes in an Effective Field Theory Extension of General Relativity, Phys. Rev. Lett. 121, 251105 (2018), arXiv:1808.08962 [gr-qc] .
  * Cano  _et al._ [2020] P. A. Cano, K. Fransen, and T. Hertog, Ringing of rotating black holes in higher-derivative gravity, Phys. Rev. D 102, 044047 (2020), arXiv:2005.03671 [gr-qc] .
  * Cano  _et al._ [2022] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, Gravitational ringing of rotating black holes in higher-derivative gravity, Phys. Rev. D 105, 024064 (2022), arXiv:2110.11378 [gr-qc] .
  * Cano  _et al._ [2023a] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, Universal Teukolsky equations and black hole perturbations in higher-derivative gravity, Phys. Rev. D 108, 024040 (2023a), arXiv:2304.02663 [gr-qc] .
  * Cano  _et al._ [2023b] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, Quasinormal modes of rotating black holes in higher-derivative gravity, (2023b), arXiv:2307.07431 [gr-qc] .
  * Weinberg [2008] S. Weinberg, Effective Field Theory for Inflation, Phys. Rev. D 77, 123541 (2008), arXiv:0804.4291 [hep-th] .
  * Mamani  _et al._ [2022] L. A. H. Mamani, A. D. D. Masa, L. T. Sanches, and V. T. Zanchin, Revisiting the quasinormal modes of the Schwarzschild black hole: Numerical analysis, Eur. Phys. J. C 82, 897 (2022), arXiv:2206.03512 [gr-qc] .
  * Hod [2008] S. Hod, Slow relaxation of rapidly rotating black holes, Phys. Rev. D 78, 084035 (2008), arXiv:0811.3806 [gr-qc] .
  * Sotiriou and Zhou [2014] T. P. Sotiriou and S.-Y. Zhou, Black hole hair in generalized scalar-tensor gravity, Phys. Rev. Lett. 112, 251102 (2014), arXiv:1312.3622 [gr-qc] .
  * Silva  _et al._ [2018] H. O. Silva, J. Sakstein, L. Gualtieri, T. P. Sotiriou, and E. Berti, Spontaneous scalarization of black holes and compact stars from a Gauss-Bonnet coupling, Phys. Rev. Lett. 120, 131104 (2018), arXiv:1711.02080 [gr-qc] .
  * Dima  _et al._ [2020] A. Dima, E. Barausse, N. Franchini, and T. P. Sotiriou, Spin-induced black hole spontaneous scalarization, Phys. Rev. Lett. 125, 231101 (2020), arXiv:2006.03095 [gr-qc] .
  * Herdeiro  _et al._ [2021] C. A. R. Herdeiro, E. Radu, H. O. Silva, T. P. Sotiriou, and N. Yunes, Spin-induced scalarized black holes, Phys. Rev. Lett. 126, 011103 (2021), arXiv:2009.03904 [gr-qc] .
  * Chandrasekhar [1985] S. Chandrasekhar,  _The mathematical theory of black holes_ (1985).
  * Davies and Reall [2022] I. Davies and H. S. Reall, Well-posed formulation of Einstein-Maxwell effective field theory, Phys. Rev. D 106, 104019 (2022), arXiv:2112.05603 [gr-qc] .
  * Drummond and Hathrell [1980] I. T. Drummond and S. J. Hathrell, QED Vacuum Polarization in a Background Gravitational Field and Its Effect on the Velocity of Photons, Phys. Rev. D 22, 343 (1980).
  * Daniels and Shore [1996] R. D. Daniels and G. M. Shore, ’Faster than light’ photons and rotating black holes, Phys. Lett. B 367, 75 (1996), arXiv:gr-qc/9508048 .
  * Dalvit  _et al._ [2001] D. A. R. Dalvit, F. D. Mazzitelli, and C. Molina-Paris, One loop graviton corrections to Maxwell’s equations, Phys. Rev. D 63, 084023 (2001), arXiv:hep-th/0010229 .
  * Chen  _et al._ [2017] S. Chen, S. Wang, Y. Huang, J. Jing, and S. Wang, Strong gravitational lensing for the photons coupled to a Weyl tensor in a Kerr black hole spacetime, Phys. Rev. D 95, 104017 (2017), arXiv:1611.08783 [gr-qc] .
  * Johnson  _et al._ [2023] J. P. Johnson, S. Jana, and S. Shankaranarayanan, A test of Einstein’s equivalence principle in ngEHT observations, (2023), arXiv:2303.13271 [gr-qc] .
  * Balakin and Zayats [2018] A. B. Balakin and A. E. Zayats, Non-minimal Einstein–Maxwell theory: the Fresnel equation and the Petrov classification of a trace-free susceptibility tensor, Class. Quant. Grav. 35, 065006 (2018), arXiv:1710.08013 [gr-qc] .
  * Chen and Jing [2013] S. Chen and J. Jing, Dynamical evolution of the electromagnetic perturbation with Weyl corrections, Phys. Rev. D 88, 064058 (2013), arXiv:1307.7459 [gr-qc] .
  * Wald [1984] R. M. Wald,  _General Relativity_ (Chicago Univ. Pr., Chicago, USA, 1984).
  * Geroch  _et al._ [1973] R. P. Geroch, A. Held, and R. Penrose, A space-time calculus based on pairs of null directions, J. Math. Phys. 14, 874 (1973).
  * Newman and Penrose [1962] E. Newman and R. Penrose, An Approach to gravitational radiation by a method of spin coefficients, J. Math. Phys. 3, 566 (1962).
  * Penrose and Rindler [2011] R. Penrose and W. Rindler,  _Spinors and Space-Time_, Cambridge Monographs on Mathematical Physics (Cambridge Univ. Press, Cambridge, UK, 2011).
  * Frolov and Novikov [1998] V. P. Frolov and I. D. Novikov, eds.,  _Black hole physics: Basic concepts and new developments_ (1998).
  * Kinnersley [1969] W. Kinnersley, Type D Vacuum Metrics, J. Math. Phys. 10, 1195 (1969).
  * Aksteiner and Andersson [2011] S. Aksteiner and L. Andersson, Linearized gravity and gauge conditions, Class. Quant. Grav. 28, 065001 (2011), arXiv:1009.5647 [gr-qc] .
  * Borissov and Fiziev [2010] R. S. Borissov and P. P. Fiziev, Exact Solutions of Teukolsky Master Equation with Continuous Spectrum, Bulg. J. Phys. 37, 065 (2010), arXiv:0903.3617 [gr-qc] .
  * Newman and Penrose [1966] E. T. Newman and R. Penrose, Note on the Bondi-Metzner-Sachs group, J. Math. Phys. 7, 863 (1966).
  * Fiziev [2010] P. P. Fiziev, Classes of Exact Solutions to the Teukolsky Master Equation, Class. Quant. Grav. 27, 135001 (2010), arXiv:0908.4234 [gr-qc] .
  * Ansorg and Panosso Macedo [2016] M. Ansorg and R. Panosso Macedo, Spectral decomposition of black-hole perturbations on hyperboloidal slices, Phys. Rev. D 93, 124016 (2016), arXiv:1604.02261 [gr-qc] .
  * Panosso Macedo [2020] R. Panosso Macedo, Hyperboloidal framework for the Kerr spacetime, Class. Quant. Grav. 37, 065019 (2020), arXiv:1910.13452 [gr-qc] .
  * Ripley [2022] J. L. Ripley, Computing the quasinormal modes and eigenfunctions for the Teukolsky equation using horizon penetrating, hyperboloidally compactified coordinates, Class. Quant. Grav. 39, 145009 (2022), arXiv:2202.03837 [gr-qc] .
  * Wald [1978] R. M. Wald, Construction of Solutions of Gravitational, Electromagnetic, Or Other Perturbation Equations from Solutions of Decoupled Equations, Phys. Rev. Lett. 41, 203 (1978).
  * Reall [2021] H. S. Reall, Causality in gravitational theories with second order equations of motion, Phys. Rev. D 103, 084027 (2021), arXiv:2101.11623 [gr-qc] .
  * Motl and Neitzke [2003] L. Motl and A. Neitzke, Asymptotic black hole quasinormal frequencies, Adv. Theor. Math. Phys. 7, 307 (2003), arXiv:hep-th/0301173 .
  * Will [2014] C. M. Will, The Confrontation between General Relativity and Experiment, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .

  *[BH]: black hole
  *[GR]: general relativity
  *[BHs]: black hole
  *[EOM]: equations of motion
  *[EOMs]: equations of motion
  *[QNM]: quasinormal mode
  *[QNMs]: quasinormal mode
  *[GWs]: Gravitational Wave
  *[EFT]: effective field theory
  *[EFTs]: effective field theory
  *[ODEs]: ordinary differential equation
  *[GB]: Gauss-Bonnet
  *[EsGB]: Einstein-scalar-Gauss-Bonnet
  *[RHS]: right hand side
  *[KG]: Klein-Gordon
  *[GW]: Gravitational Wave
  *[BL]: Boyer-Lindquist
  *[GHP]: Geroch-Held-Penrose
  *[NP]: Newman-Penrose
  *[PDE]: partial differential equation
  *[ODE]: ordinary differential equation
  *[BCs]: boundary condition
  *[BVP]: boundary value problem
  *[LHS]: left hand side
  *[PDEs]: partial differential equation
  *[BVPs]: boundary value problem