# EFT corrections to scalar and vector quasinormal modes of rapidly rotating
black holes

Filipe S. Miguel<EMAIL_ADDRESS>Department of Applied Mathematics and
Theoretical Physics, University of Cambridge,
Wilberforce Road, Cambridge CB3 0WA, United Kingdom

###### Abstract

Quasinormal modes characterize the final stage of a black hole merger. In this
regime, spacetime curvature is high, these modes can be used to probe
potential corrections to general relativity. In this paper, we utilize the
effective field theory framework to compute the leading order correction to
massless scalar and electromagnetic quasinormal modes. Proceeding
perturbatively in the size of the effective field theory length scale, we
describe a general method to compute the frequencies for Kerr black holes of
any spin. In the electromagnetic case, we study both parity even and parity
odd effective field theory corrections, and, surprisingly, prove that the two
have the same spectrum. Furthermore, we find that, the corrected frequencies
separate into two families, corresponding to the two polarizations of light.
The corrections pertaining to each family are equal and opposite. Our results
are validated through several consistency checks.

BC

    boundary condition
BH

    black hole
BL

    Boyer-Lindquist
BVP

    boundary value problem
EsGB

    Einstein-scalar-Gauss-Bonnet
EFT

    effective field theory
EH

    Einstein-Hilbert
EM

    Einstein-Maxwell
EOM

    equations of motion
GB

    Gauss-Bonnet
GHP

    Geroch-Held-Penrose
GR

    general relativity
GW

    Gravitational Wave
KG

    Klein-Gordon
LHS

    left hand side
NP

    Newman-Penrose
ODE

    ordinary differential equation
PDE

    partial differential equation
QNM

     quasinormal mode
RHS

    right hand side

###### Contents

  1. I Introduction
  2. II Background Material
    1. II.1 The Kerr spacetime
    2. II.2 The GHP formalism
    3. II.3 Master wave equation
    4. II.4 Quasinormal Modes
    5. II.5 Parity transformations
    6. II.6 Hodge duals and p forms
    7. II.7 The vacuum Maxwell equations
  3. III Corrections to scalar modes
    1. III.1 Parity even equations of motion
    2. III.2 On parity odd EFT corrections
    3. III.3 Results
  4. IV Corrections to electromagnetic modes
    1. IV.1 Equations of motion
    2. IV.2 Parity even corrections
      1. IV.2.1 Eikonal limit
      2. IV.2.2 Schwarzschild limit
      3. IV.2.3 Generalization to Kerr
    3. IV.3 Mixed parity corrections
      1. IV.3.1 Schwarzschild limit
      2. IV.3.2 Generalization to Kerr
    4. IV.4 Results
      1. IV.4.1 Schwarzschild limit
      2. IV.4.2 Kerr results
  5. V Discussion
  6. VI Acknowledgements
  7. A Numerical methods
    1. A.1 Pseudo-spectral methods
    2. A.2 Solving the background equation
    3. A.3 Solving the perturbation equation

##  I Introduction

Studying black hole (BH) spacetimes has proven to be a very fruitful endeavor
in understanding the merits and limitations of Einstein’s theory of general
relativity (GR). These powerful astrophysical objects cause spacetime to
heavily curve around them, and thus, we expect possible deviations from GR to
be noticeable there. The no-hair theorems imply that the mathematical
description of these objects is remarkably simple. In a vacuum, asymptotically
flat spacetime, all stationary BH solutions are described by the two-parameter
Kerr family of metrics. The mass and angular momentum $\left(M,a\right)$ of
the BH are sufficient to specify these highly symmetric solutions of the
Einstein equations. In the limit $a\rightarrow 0$, we obtain the only static,
spherically symmetric solution of the equations of motions, known as the
Schwarzschild solution. However, real-world scenarios are dynamic, and thus,
solutions will necessarily differ from the stationary ones discussed above.

The observation of Gravitational Waves by the LIGO/VIRGO experiment offers a
key opportunity to experimentally test the dynamical evolution of BH
spacetimes. In particular, by observing the collision between two BHs, we can
understand how generic BH solutions evolve into the stationary Kerr spacetime.
The collision of two BHs is characterized by three phases. Initially, during
the inspiral phase, the BHs orbit each other, with motion that can be
described using Newtonian physics. Then, during the merger phase, the BHs
effectively collide and form a single BH. This phase is challenging to study
analytically, and a numerical relativity approach must be adopted. Finally,
during the BH ringdown, the final BH decays into stationarity, ’vibrating’
around a Kerr solution. This phase can be modelled by understanding linear
perturbations of the EOM.

The study of the linear perturbations of Schwarzschild was pioneered in [1],
where the authors derived the Regge-Wheeler equation, which describes odd
parity (axial) perturbations of a Schwarzschild spacetime. This work was
complemented in [2] with the Zerilli equation, describing even parity (polar)
perturbations of Schwarzschild. Concurrently, Teukolsky demonstrated that in a
Kerr background, most fields with linear EOMs are described by two coupled
wave-like ordinary differential equations [3]. The Schwarzschild limit of this
is known as the Bardeen-Press equation [4], which does not coincide with the
Zerilli or Regge-Wheeler equation. It was later proved that all three
equations are related to each other by ordinary and generalized Darboux
transformations [5, 6]. This implies that the polar and axial EOM have the
same spectrum, a conclusion that does not necessarily generalize to all
physical systems, as we will discuss later. Therefore, conveniently,
characterising solutions to the Teukolsky equation is sufficient to understand
the behaviour of most linear fields in a Kerr/Schwarzschild background.

Solving these equations in general is no simple task, however, similar to
quantum mechanical systems, we can decompose general solutions into a set of
normal modes and frequencies - the eigenvector/eigenvalue pairs of the system.
BHs are, by definition, absorbent objects: the energy of a field in a BH
background is gradually dissipated into its interior. Vishveshwara [7] was the
first to understand that this implies the frequencies must be complex, with
the imaginary part controlling the damping of the fields. Thus, these are
quasinormal modes as was coined in [8]. See [9, 10, 11, 12] for a review.

Obtaining the values of the QNM frequencies is a non-trivial problem.
Analytically, we can use a WKB approach to approximate the values of large QNM
frequencies. In fact, it can be demonstrated that the frequency and damping of
these QNMs are related to the frequency and Lyapunov exponent of circular null
geodesics around the BH, as seen in [13, 14, 15]. This connection will be
vital to make progress in section IV.2.1. However, to achieve accurate values,
we must resort to a numerical approach. Following a strategy to compute the
spectrum of the molecular hydrogen ion, Leaver devised the first method that
yields accurate results across most of the Kerr parameter space [16]. However,
the approach loses accuracy in the extremal limit ($a\rightarrow M$). Modern
approaches project the EOM onto a discrete grid, finding the QNM spectrum with
a pseudospectral approach, as detailed in [17, 18, 19, 20, 21].

Since QNM frequencies are complex-valued, the EOM cannot be self-adjoint.
Additional complexity arises from the fact that, in general, the QNMs are not
a complete basis for the solutions of the EOM; for an arbitrary linear
perturbation, there will be a residual in the form of a polynomial tail.
Nevertheless, we expect that most of the ringdown waveform can be described
with QNM contributions, see [22] for details. By fitting the data to a finite
sum of the slowest decaying modes, we can obtain the values of the first few
QNM frequencies. This is known in the literature as BH spectroscopy. In a Kerr
background, the frequencies are fully specified by the mass and spin of the
BH. Theoretically, by extracting the real and imaginary parts of a single QNM
with sufficient precision, we have enough data to fully specify the parameters
of the BH. Conversely, by detecting two or more modes, we can test GR as a
theory [23]. If the detected spectrum is inconsistent with the predictions
from GR, we need to consider corrections to GR. Modern approaches use the GWs
from the inspiral and merger phases to determine the BH parameters, and then
check their consistency with the ringdown spectrum [24, 25, 26].

To validate GR as a theory, we need a concrete framework to parametrize
possible deviations. If we knew the full form of the true theory of quantum
gravity, we could make physical predictions at low energy scales by
integrating out the high energy degrees of freedom. This process maps the
Lagrangian into an infinite sum of low-energy interactions. Using dimensional
analysis, we know that each spacetime derivative must be preceded by some
length scale $L$. Thus, we can organize the sum in ascending order of powers
of $L$, its scaling dimension, see [27, 28, 29]. Truncating the sum at some
order we obtain an effective field theory (EFT) description of the system. The
approximation is valid as long as $L$ is the smallest length scale of the
theory. To consider values of $L$ closer to the spacetime length scale, we
must include more terms in the sum. [30] provides compelling evidence that at
a classical level, the difference between the solutions to EFT and UV EOM is
dominated by the scaling dimension of the next term in the series. Remarkably,
making simple assumptions about the theory (e.g., gauge invariance and
locality), the form of the interactions becomes heavily constrained. Thus, up
to a given accuracy, the space of EFT is fully parametrized by a set of
finitely many coupling constants. This suggests a bottom-up approach. We
include in our Lagrangian all the physically reasonable terms up to a given
order, extract physical predictions, and use experiments to constrain the
coupling constant values. Then, we eliminate all UV theories that lead to
coupling parameters that violate these constraints. Thus, even if we do not
observe any deviations from GR, we can constrain the UV theory space.

Calculating the impact of gravitational EFT corrections on QNMs is tricky.
Take some stationary BH solution of the uncorrected GR EOM. QNMs are encoded
in dynamical perturbations of this background, i.e. GWs. Gravitational EFTs
will correct both the stationary background metric and the GWs around it.
Crucially, the correction to the GWs will have the same order of magnitude as
the correction to the stationary background. Thus, the EFT corrections to the
QNMs EOM depend on the EFT correction to the background metric and pure
corrections to the linearised GR equations. In the Schwarzschild limit, [31]
computed the correction to the background spacetime, and corresponding
correction to gravitational QNMs. More recently, Cano et al. studied the
corrections to Kerr QNMs. In [32], they obtained the correction to the
background metric for BHs with $a\lesssim 0.7M$, and used the result to obtain
the shifts in the QNMs of a massless scalar field. Then, in [33, 34, 35], the
same group obtained the corrections to gravitational QNMs of not rapidly
rotating Kerr BHs, under parity even and parity odd EFT terms. Crucially, all
of the results are perturbative in $a/M$ (where $a$ is the BH spin parameter),
with the method breaking down for rapidly rotating BHs. In fact, in [34], the
authors argue that the method breaks down when considering BHs with $a\gtrsim
0.7M$. Astrophysical BHs can have $a/M$ much closer to $1$, thus, it is
essential to develop an approach that is non-perturbative in $a/M$. This is
fundamentally more difficult, because the EOM cannot be separated into radial
and angular ODEs, the QNMs are eigenvalues of two dimensional partial
differential equations. Throughout this paper, we develop a framework to
perform this calculation by studying similar, albeit simpler, theories. This
is in preparation for future work that will generalize the method to the
gravitational case. We will study leading order EFT corrections to scalar and
electromagnetic QNMs of Kerr BHs of any spin.

We begin by studying the scalar QNMs of an EFT that couples gravity with a
massless scalar field $\Phi$. We require that $\Phi$ is small and that, in the
$\Phi\rightarrow 0$ limit, we recover classical GR solutions. This implies
that, to leading order, we must restrict to EFT corrections that lead to EOM
that are linear in $\Phi$. Additionally, for simplicity, we restrict to terms
that are parity invariant.111To achieve full generality, we must include
parity odd corrections, however, we will argue in section III.2 that our
method can be swiftly adapted to study this case. After suitable field
redefinitions, the leading order EFT can be written as follows [36]:

$S=\frac{1}{16\pi}\int_{\mathcal{M}}\differential^{4}x\sqrt{-g}\left[R-\frac{1}{2}\nabla_{a}\Phi\nabla^{a}\Phi\,+L^{2}f(\Phi)\,\mathcal{G}\right]\,,$
(1)

where

$\mathcal{G}:=R_{abcd}R^{abcd}-4R_{ab}R^{ab}+R^{2}\,,$ (2)

is the Gauss-Bonnet (GB) term, and $L$ is the EFT length scale. This is
commonly referred to, in the literature, as the Einstein-scalar-Gauss-Bonnet
(EsGB) theory. By varying the action with respect to $\Phi$, we derive the
EOM:

$\Box\Phi=-L^{2}f^{\prime}(\Phi)\mathcal{G}\,.$ (3)

(a)

(b)

(c)

(d)

Figure 1: This image summarizes the EFT corrections to low $\ell$ QNMs of a
massless scalar field. In panels a) and b) we have the real and imaginary
parts of the background QNM frequencies, whereas panels c) and d) show the EFT
correction. To generate these plots, we computed the QNMs for BHs with
$a\lesssim 0.999M$ in most cases, or $a\lesssim 0.999975M$ for $\ell=m\neq 0$
modes. We found that, in general, the imaginary part of the background
frequency increases with $\ell$. Thus, the lowest lying QNMs should have
$\ell\rightarrow\infty$. This is consistent with the literature, (see table 1
of [37]), but is contrary to the tendency for gravitational QNMs. As expected,
in the Schwarzschild limit, both the background and correction frequencies,
are independent of $m$, with the degeneracy being broken for $a>0$. In the
near extremal limit, for $\ell=m\neq 0$ modes $M\omega_{0}\rightarrow\ell/2$,
as predicted in [38]. This feature extends to the EFT correction:
$\omega_{1}\rightarrow 0$ as $a\rightarrow M$.

(a)

(b)

(c)

(d)

Figure 2: Here, we summarize the EFT corrections to the lowest lying
electromagnetic QNMs. The real and imaginary parts of the background QNM
frequencies can be seen in panels a) and b) whereas the correction can be
found in panels c) and d). As expected, in the Schwarzschild limit, modes are
independent of $\ell$. Contrary to the scalar case, for most of the parameter
space, $\mathcal{I}m(\omega_{0})$ decreases with $\ell$, making the mode with
$\ell=1$ the lowest lying mode. This tendency is broken at $a\approx 0.96M$,
where the $\ell=2$ modes become the lowest lying QNM. As in the scalar case,
for $\ell=m$ modes, in the limit $a\rightarrow M$, we have
$M\omega_{0}\rightarrow\ell/2$. This is in agreement with the prediction in
[38]. Furthermore, as seen in the scalar case, the corresponding correction
obeys $M\omega_{1}\rightarrow 0$.

Different choices of $f(\Phi)$ have different consequences. If $f(\Phi)$ is
constant, the right hand side (RHS) vanishes, and we recover the massless
Klein-Gordon (KG) equation. Conversely, if $f(\Phi)\sim\Phi$, the standard GR
BH solutions no longer solve the gravitational EOM, resulting in BHs solutions
with scalar hair [39]. To avoid this, we must consider second order:
$f(\Phi)\sim\Phi^{2}$. For this theory, $\Phi=0$ solves the scalar EOM, and
thus, GR BHs are in the solution space. Take $f(\Phi)=\lambda\Phi^{2}/2$ for
some $\mathcal{O}(1)$ dimensionless coupling parameter $\lambda$, the EOM for
the scalar field becomes:

$\Box\Phi+L^{2}\,\lambda\,\mathcal{G}\,\Phi=0\,.$ (4)

This equation is analogous to the massive KG equation
$\left(\Box-m^{2}\right)\Phi=0$. In a flat background, when $m^{2}>0$, $\Phi$
is stable. However, if $m^{2}<0$, the KG field is susceptible to the well-
known Tachyonic instability. The GB coupling acts as a spacetime varying
effective mass term, thus, we expect stability to be governed by its sign.
[40] examined equation (4) in a Schwarzschild background. The authors
discovered that for $\lambda<0$, the effective mass is positive, and no hairy
BH solutions can form; all static solutions default to the standard
Schwarzschild family of BHs. On the other hand, for $\lambda>0$, BHs with
sufficiently small mass ($M$) can be unstable. Indeed, the authors show that
for certain mass ranges obeying $\lambda\left(\frac{L}{M}\right)^{2}\gtrsim
0.6$, classical BH solutions dynamically evolve into hairy BHs in a process
known as BH scalarization. When extending this to BHs with arbitrary spin, we
get surprising results. [41, 42] showed that if the BH is sufficiently close
to extremality, $\mathcal{G}$ can be negative, reversing the sign of the
effective mass. In this case, the Tachyonic instability emerges for
$\lambda<0$. Again, this only occurs for large enough
$\absolutevalue{\lambda}(L/M)^{2}$. All of these results occur outside the
regime of validity of the EFT, and thus, could be suppressed or mitigated by
higher-order derivative couplings. Crucially, using the action (1), we are
only probing corrections to observables that are at most quadratic in $L/M$,
or equivalently, linear in $\lambda$. All these papers investigate the
instability by time-evolving equation (4), and inferring an instability time
scale. As far as we know, there are no results in the literature following a
standard QNM approach.

In section III, we will calculate the corrections to scalar QNMs for BHs of
any spin, working perturbatively in $L/M$. The frequencies take the form:

$\omega=\omega_{0}\,+\,\lambda\left(\frac{L}{M}\right)^{2}\,\omega_{1}\,.$ (5)

For the reader’s convenience, we have summarized the main results in figure 1.
We will discuss the results in detail in section III.3. For now, we emphasize
that these plots reveal non-trivial behaviour for $a/M>0.9$ that would be
difficult to capture via an approach perturbative in $a/M$.

The electromagnetic sector bears a closer resemblance to the physically
relevant gravitational case. In both cases, there are only two degrees of
freedom, encoded by two electromagnetic / GW polarizations. In the GR limit,
the two polarizations have the same QNM spectrum, i.e. they are isospectral
[43]. Thus, by understanding the corrections to electromagnetic QNMs, we can
develop intuition that will be vital in characterizing the gravitational
sector in the future. Consider the following, EFT action:

$S=\frac{1}{16\pi}\int_{\mathcal{M}}\differential^{4}x\sqrt{-g}\,\left[R-\frac{1}{4}\left(F_{ab}F^{ab}-\frac{1}{3}L^{2}\left(\lambda_{(e)}\,\mathcal{C}_{abcd}F^{ab}F^{cd}+\lambda_{(o)}\,\mathcal{C}_{abcd}\left(\star
F\right)^{ab}F^{cd}\right)\right)\right]\,.$ (6)

Here, $L$ is the EFT length scale, $F$ represents the electromagnetic tensor,
$\mathcal{C}$ denotes the Weyl tensor, $\star$ denotes the hodge dual
operation (see section II.6), and $\lambda_{(e)}/\lambda_{(o)}$ are
dimensionless coupling constants. The terms proportional to $\lambda_{(e)}$
and $\lambda_{(o)}$ are, respectively, even / odd under parity
transformations, see section II.5 for details in this. Working perturbatively
in the size of $F$, and the EFT length scale, we can show that, up to field
redefinitions, this is the leading order of the most general EFT that couples
electromagnetism and gravity, see e.g. [44]. The action includes terms up to
$\mathcal{O}(L/M)^{2}$, thus, as in the scalar case, our results are only
valid to that order. Higher order contributions are outside the regime of
validity of the EFT.

Setting $\lambda_{(o)}=0$, the action is invariant under parity
transformations. This particular regime has been extensively studied in the
literature. [45, 46] focused on the propagation of light in Schwarzschild and
Kerr backgrounds, respectively. They concluded that the propagation properties
are tied to the photon polarisation, rendering spacetime into a birefringent
medium. The observational implications of this have been studied in [46, 47,
48, 49]. [50] generalized this result to various algebraically special
spacetimes. Notably, for Petrov type-D spacetimes, the authors found two
effective metrics that characterize the propagation of each polarization of
light. Because QNMs are related to light propagation in the near extremal
limit, we expect that the QNM spectrum will no longer be polarization
independent. [51] confirmed this was the case on a Schwarzschild background by
explicitly computing the QNMs across a wide range of $\lambda$ values.

When considering parity invariant EFT corrections, we can find a basis for
electromagnetic modes, where they are eigenstates of the parity operator.
Thus, it is natural, to associate these with the two possible polarizations.
We get the parity even / polar polarization and the parity odd / axial
polarization. However, when we consider parity odd EFT corrections, the parity
invariance of the action is broken, and thus, parity even modes will be mixed
with parity odd modes, see [33]. In section IV, we will consider, both, parity
even and parity odd EFT corrections. Working perturbatively in $L/M$, we will
prove that the shift in the QNM spectrum due to parity odd corrections is
exactly the same as the one for parity even ones. Furthermore, we find that
each background QNM frequency has two possible corrections, related by a minus
sign.

$\omega^{\pm}=\omega_{0}\pm(-1)^{\ell+m+1}\left(\frac{L}{M}\right)^{2}\absolutevalue{\lambda}\,\omega_{1}\,+\mathcal{O}\left(\frac{L}{M}\right)^{3}\,,$
(7)

where:

$\absolutevalue{\lambda}=\sqrt{\lambda_{(e)}^{2}+\lambda_{(o)}^{2}}\,.$ (8)

Cano et. al, obtained similar results for the gravitational case [33, 34, 35].
The authors found similar mixing between parity even and parity odd
corrections. Possibly because the EFT corrects the background spacetime, the
authors found that $\omega_{1}^{+}\neq-\omega_{1}^{-}$. The corrected QNM
frequencies are not centred around $\omega_{0}$. For the reader’s convenience,
we summarized the main results in figure 2. As before, the plots reveal non-
trivial behaviour for $a/M\gtrsim 0.75$ that would be difficult to capture via
a perturbative approach. We will discuss these results in detail in section
IV.4.2.

In the context of QED in a curved background, when we integrate out the
electron degrees of freedom, we get the parity even part of the action (6). In
this context, the value of $\lambda_{(e)}$ is known [45]:

$\frac{1}{3}\lambda_{(e)}L^{2}=\frac{1}{1440\pi^{2}}\left(\frac{M_{p}^{2}}{m_{e}M_{\odot}}\right)^{2}M_{\odot}^{2}\,\sim\,4.8\cdot
10^{-36}M_{\odot}^{2}\,.$ (9)

Here, $M_{p}$ is the Planck mass, $m_{e}$ is the mass of an electron, and
$M_{\odot}$ is the mass of the Sun. For a solar mass BH, we anticipate the
effects to be extremely small, with corrections only becoming relevant for BHs
of very small masses. However, there might be contributions to $\lambda_{(e)}$
that are unrelated to QED, thus its true value could be much larger. Despite
that, this action is still interesting as a toy model for the pure gravity
case.

In the calculations of this paper, we followed the $(-,+,+,+)$ signature
convention. We perform our calculations using geometric units, i.e.
$(G_{N}=c=1)$. This implies that dimensions of any quantity can be expressed
as the power of some length scale. Furthermore, we use Latin indices
$a,b,c,\cdots$ for abstract index notation, Greek indices
$\mu,\nu,\sigma,\cdots$ for quantities in a coordinate basis, and
parenthesized Greek indices $(\alpha),(\beta),(\gamma),\cdots$ for quantities
in a null tetrad basis.

##  II Background Material

###  II.1 The Kerr spacetime

Throughout this paper, we will work on a Kerr background. This is an
asymptotically flat solution of the vacuum Einstein equations, that describes
a rotating BH. It is fully specified by two parameters, $(M,a)$, representing
the mass and angular momentum of the BH, respectively. We will work with
modified Boyer-Lindquist (BL) coordinates $(t,r,x,\phi)$, that differ from
standard BL coordinates by the definition $x=\cos\theta$. The line element
reads:

$\differential s^{2}=-\differential t^{2}\,-\,\frac{2Mr\left(\differential
t-a\,\Delta_{x}\,\differential\phi\right)^{2}}{\Sigma}\\\
+\frac{\Sigma\,\differential
r^{2}}{\Delta_{r}}\,+\,\frac{\Sigma\,\differential
x^{2}}{\Delta_{x}}\,+\,\Delta_{x}\left(r^{2}+a^{2}\right)\differential\phi^{2},$
(10)

where:

$\displaystyle\Sigma=r^{2}+a^{2}x^{2},$ (11)
$\displaystyle\Delta_{r}=r^{2}-2Mr+a^{2}=(r-r_{+})(r-r_{-}),$
$\displaystyle\Delta_{x}=1-x^{2},$ $\displaystyle
r_{\pm}=M\pm\sqrt{M^{2}-a^{2}},$ $\displaystyle t\in\mathbb{R},\quad
r\in(r_{+},+\infty),\quad x\in(-1,1),\quad\phi\in(0,2\pi).$

We will restrict our study to BHs such that $M\in\mathbb{R}^{+}$. In the limit
$a\rightarrow 0$, the spacetime is spherically symmetric, and the line element
reduces to the Schwarzschild metric. On the other hand, at $a=M$, we have the
extremal limit. If we allowed $a>M$, $\Delta_{r}$ would have no real roots,
and the metric would describe a naked singularity. Therefore, we restrict to
$0<a<M$. The causal structure of the BH exterior can be seen in figure 3.
$\mathcal{H}^{+}$ and $\mathcal{H}^{-}$ are the BH and white hole horizons,
respectively, whereas $\mathcal{J}^{\pm}$ represents future and past null
infinity. Using Kerr coordinates, the spacetime can be analytically extended
to the BH interior (see [52]). At $r=r_{-}$ there is a Cauchy horizon
($\mathcal{CH}$). The surface gravity $\kappa$ of $\mathcal{H}^{+}$ is:

$\kappa=\frac{1}{2\,a}\Delta_{r}^{\prime}(r_{+})\,\Omega_{+}=\frac{r_{+}-r_{-}}{2\left(a^{2}+r_{+}^{2}\right)}\,,$
(12)

where:

$\Omega_{+}=\frac{a}{a^{2}+r_{+}^{2}}\,.$ (13)

This parameter describes the angular velocity of $\mathcal{H}^{+}$ as measured
by an observer at infinity. In the extremal limit, $r_{-}\rightarrow r_{+}$,
and $\kappa\rightarrow 0$. In this limit, we no longer have a Cauchy horizon,
and $\mathcal{H}$ becomes a doubly degenerate Killing horizon. This behavior
makes calculations in the extremal limit increasingly difficult.

Finally, we define the dimensionless angular momentum:

$J=\frac{a}{M}$ (14)

If we replace $a$ by $J$ in the metric, a change $M\rightarrow\alpha M$ is
equivalent to the coordinate transformation $r\rightarrow
r/\alpha\,,t\rightarrow t/\alpha$. Consequently, all dimensionless background
quantities must be independent of $M$. $J$ is the only physically relevant
parameter.

$i^{+}$$i^{-}$$i^{0}$$\cal{J}^{+}$$r\rightarrow+\infty\quad
t\rightarrow+\infty$$\mathcal{H}^{+}$$r\rightarrow r_{+}\quad
t\rightarrow+\infty$$\mathcal{H}^{-}$$r\rightarrow r_{+}\quad
t\rightarrow-\infty$$\cal{J}^{-}$$r\rightarrow+\infty\quad
t\rightarrow-\infty$

Figure 3: Penrose diagram for the exterior of a Kerr BH.

###  II.2 The GHP formalism

The Geroch-Held-Penrose (GHP) formalism [53] is a more covariant version of
the Newman-Penrose (NP) approach [54, 55]. This formalism is useful to derive
compact EOM for various fields, especially when working within the context of
algebraically special spacetimes.

As in the NP formalism, the core idea is to project all relevant quantities
and equations into a complex null tetrad basis $e^{a}{}_{(\alpha)}$. To define
this structure, we begin by selecting two null curve congruences, $\gamma_{l}$
and $\gamma_{n}$, with tangents $l^{a}$ and $n^{a}$ satisfying $l\cdot n=-1$.
Subsequently, we choose two spacelike vector fields $X^{a}$ and $Y^{a}$
orthogonal to both $l^{a}$ and $n^{a}$, such that:

$\displaystyle X_{a}X^{a}=Y_{a}Y^{a}=1/2\,,$ (15) $\displaystyle
X_{a}Y^{a}=0\,.$

These are then used to define complex vector fields $m^{a}$ and $\bar{m}^{a}$:

$\displaystyle m^{a}$ $\displaystyle:=X^{a}+iY^{a}\,,$ (16)
$\displaystyle\bar{m}^{a}$ $\displaystyle:=X^{a}-iY^{a}\,.$

Throughout this paper, the bar $(\,\bar{}\,)$ operation will denote complex
conjugation. As a direct consequence of equation (15), we find
$m^{2}=\bar{m}^{2}=0$ and $m\cdot\bar{m}=1$. Finally, the null tetrad will be
$e^{a}{}_{(\alpha)}:=(l^{a},n^{a},m^{a},\bar{m}^{a})_{(\alpha)}$. The
normalization conditions reduce to:

$g_{ab}\,e^{a}{}_{(\alpha)}\,e^{b}{}_{(\beta)}\,=\,\hat{\eta}_{(\alpha)(\beta)}\,,$
(17)

with:

$\hat{\eta}_{(\alpha)(\beta)}=\begin{pmatrix}0&-1&0&0\\\ -1&0&0&0\\\
0&0&0&1\\\ 0&0&1&0\end{pmatrix}_{(\alpha)(\beta)}\,.$ (18)

After fixing the directions of $\gamma_{l}$ and $\gamma_{n}$, there are still
multiple consistent choices of $e_{(\alpha)}$, related to each other by a
continuous transformation. We understand this redundancy as the gauge freedom
of the formalism. On one hand, we may reparametrize $\gamma_{l}$ and
$\gamma_{n}$, such that:

$l^{a}\rightarrow r\,l^{a},\quad n^{a}\rightarrow\frac{1}{r}\,n^{a}\,,$ (19)

where $r>0$ may vary across the spacetime. On the other hand, we can rotate
the spacelike basis used to define $m^{a}$ by a real angle $\theta$, which,
again, need not be a constant:

$m^{a}\rightarrow e^{i\theta}\,m^{a},\quad\bar{m}^{a}\rightarrow
e^{-i\theta}\,\bar{m}^{a}\,.$ (20)

The two transformations can be encoded in a single scalar field
$\lambda=\sqrt{r}e^{i\theta/2}$:

$\displaystyle l^{a}\rightarrow\lambda\bar{\lambda}\,l^{a},\quad$
$\displaystyle n^{a}\rightarrow\lambda^{-1}\bar{\lambda}^{-1}\,n^{a},$ (21)
$\displaystyle m^{a}\rightarrow\lambda\bar{\lambda}^{-1}\,m^{a},\quad$
$\displaystyle\bar{m}^{a}\rightarrow\lambda^{-1}\bar{\lambda}\,\bar{m}^{a}$

In general, if under (21), $\eta$ transforms as:

$\eta\rightarrow\lambda^{p}\bar{\lambda}^{q}\eta\,,$ (22)

we say that $\eta$ is GHP-covariant, with type (or weight) ${p,q}$. We deduce
that the tetrad components have type:

$\displaystyle l^{a}:$ $\displaystyle\\{1,1\\}\,,$ (23) $\displaystyle n^{a}:$
$\displaystyle\\{-1,-1\\}\,,$ $\displaystyle m^{a}:$
$\displaystyle\\{1,-1\\}\,,$ $\displaystyle\bar{m}^{a}:$
$\displaystyle\\{-1,1\\}\,.$

To preserve GHP covariance, we cannot add two quantities with different
weights. Conversely, multiplying a $\\{p,q\\}$ scalar by a
$\\{p^{\prime},q^{\prime}\\}$ scalar is allowed, yielding a
$\\{p+p^{\prime},q+q^{\prime}\\}$ scalar. To further simplify the notation, we
introduce the prime $(\,^{\prime}\,)$ operation. This effectively interchanges
instances of $l^{a}$ and $n^{a}$, and swaps instances of $m^{a}$ and
$\bar{m}^{a}$:

$\left(l^{a}\right)^{\prime}=n^{a},\mkern
4.0mu\left(n^{a}\right)^{\prime}=l^{a},\mkern
4.0mu\left(m^{a}\right)^{\prime}=\bar{m}^{a},\mkern
4.0mu\left(\bar{m}^{a}\right)^{\prime}=m^{a}\,.$ (24)

If $\eta$ is of type $\\{p,q\\}$, then $\eta^{\prime}$ will be of type
$\\{-p,-q\\}$. Conversely, $\bar{\eta}$ will be of type $\\{q,p\\}$. Both
operations commute, and
$\left(\eta^{\prime}\right)^{\prime}=\bar{\bar{\eta}}=\eta$. It is also useful
to define the spin and boost weights of $\eta$ as $s=\frac{1}{2}(p-q)$ and
$p=\frac{1}{2}(p+q)$, respectively. Under complex conjugation, the spin
weight’s sign is flipped, whereas the boost weight remains unchanged. Under a
$(\,^{\prime}\,)$ transformation, both quantities have their signs inverted.

To make progress, we must define the GHP equivalent of the NP spin
coefficients [54]. These can be divided into two groups222These quantities
differ from those in [53] by a $(-)$ sign. This is because we follow the
$(-,+,+,+)$ signature convention, whereas [53] use a $(+,-,-,-)$ signature.
See appendix E of [56] for further details.:

$\displaystyle\kappa:=-m^{b}l^{a}\nabla_{a}l_{b}\,,$
$\displaystyle\sigma:=-m^{b}m^{a}\nabla_{a}l_{b}\,,$ (25)
$\displaystyle\rho:=-m^{b}\bar{m}^{a}\nabla_{a}l_{b}\,,$
$\displaystyle\tau:=-m^{b}n^{a}\nabla_{a}l_{b}\,,$

and

$\displaystyle\beta:=-1/2\left(n^{b}m^{a}\nabla_{a}l_{b}-\bar{m}^{b}m^{a}\nabla_{a}m_{b}\right)\,,$
(26)
$\displaystyle\varepsilon:=-1/2\left(n^{b}l^{a}\nabla_{a}l_{b}-\bar{m}^{b}l^{a}\nabla_{a}m_{b}\right)\,,$

together with the corresponding primed versions. One can show that the scalars
in (25) are GHP-covariant, with types:

$\kappa:\\{3,1\\}\quad\sigma:\\{3,-1\\}\quad\rho:\\{1,1\\}\quad\tau:\\{1,-1\\}.$
(27)

Conversely, because $\lambda$ is spacetime varying, the scalars in (25) do not
transform according to (22), gaining an additional affine contribution:

$\displaystyle\beta$
$\displaystyle\rightarrow\lambda\bar{\lambda}^{-1}\,\beta+\bar{\lambda}^{-1}\,m^{a}\nabla_{a}\lambda\,,$
(28) $\displaystyle\varepsilon$
$\displaystyle\rightarrow\lambda\bar{\lambda}\,\varepsilon+\bar{\lambda}\,l^{a}\nabla_{a}\lambda\,.$

A similar problem arises when taking the covariant derivative of GHP-covariant
quantities. An affine contribution will stem from the derivatives of $\lambda$
and $\bar{\lambda}$. Consider
$l^{a}\nabla_{a}A_{(l)}:=l^{a}\nabla_{a}\left(l^{b}A_{b}\right)$. Under (21)
we obtain:

$l^{a}\nabla_{a}A_{(l)}\rightarrow\lambda^{2}\bar{\lambda}^{2}\,l^{a}\nabla_{a}A_{(l)}+\lambda\bar{\lambda}\,A_{(l)}\,l^{a}\nabla_{a}\left(\lambda\bar{\lambda}\right)\,.$
(29)

Comparing equations (28) and (29), we deduce that for the quantity
$\left(l^{a}\nabla_{a}-\varepsilon-\bar{\varepsilon}\right)A_{(l)}$, the
affine contributions are cancelled. Thus, this quantity is GHP-covariant, with
weight ${1+1,1+1}$. Therefore, to preserve GHP covariance under
differentiation, we need to supplement the derivative operators with
appropriate combinations of the spin coefficients. This is entirely analogous
to the introduction of Christoffel symbols in the definition of standard
covariant derivatives. We define four derivative operators:

$\displaystyle\text{\thorn}\,\eta$
$\displaystyle=\left(l^{a}\nabla_{a}-p\,\varepsilon-q\,\bar{\varepsilon}\right)\,\eta:\quad\quad$
$\displaystyle\\{1,1\\}\,,$ (30) $\displaystyle\text{\thorn}^{\prime}\,\eta$
$\displaystyle=\left(n^{a}\nabla_{a}+p\,\varepsilon^{\prime}+q\,\bar{\varepsilon}^{\prime}\right)\,\eta:\quad\quad$
$\displaystyle\\{-1,-1\\}\,,$ (31) $\displaystyle\eth\,\eta$
$\displaystyle=\left(m^{a}\nabla_{a}+p\,\beta-q\,\bar{\beta}^{\prime}\right)\,\eta:\quad\quad$
$\displaystyle\\{1,-1\\}\,,$ (32) $\displaystyle\eth^{\prime}\,\eta$
$\displaystyle=\left(\bar{m}^{a}\nabla_{a}-p\,\beta^{\prime}+q\,\bar{\beta}\right)\,\eta:\quad\quad$
$\displaystyle\\{-1,1\\}\,,$ (33)

where the numbers on the right denote the GHP weight. þ,
$\text{\thorn}^{\prime}$ respectively raise and lower the boost-weight,
leaving the spin-weight unchanged. Conversely, $\eth$ and $\eth^{\prime}$ act
as raising and lowering operators for the spin-weight, leaving the boost-
weight unchanged. Finally, we can define a new (doubly) covariant derivative
operator that transforms covariantly under both coordinate and tetrad
transformations:

$\Theta_{a}:=-\,n_{a}\,\text{\thorn}\,-\,l_{a}\text{\thorn}^{\prime}\,+\,\bar{m}_{a}\eth\,+\,m_{a}\eth^{\prime}\,.$
(34)

This operator obeys a Leibniz rule. For $\alpha,\beta$ GHP-covariant
quantities:

$\Theta_{a}\left(\alpha\beta\right)=(\Theta_{a}\alpha)\beta\,+\,\alpha(\Theta_{a}\beta)\,.$
(35)

In 4 dimensions, the algebraic symmetries of the Weyl tensor allow for 10
degrees of freedom. Projecting it into the null tetrad, we can encode them in
5 complex scalars:

$\displaystyle\psi_{0}$ $\displaystyle:=C_{1313}=\psi_{4}^{\prime}:$
$\displaystyle\\{4,0\\},$ (36) $\displaystyle\psi_{1}$
$\displaystyle:=C_{1213}=\psi_{3}^{\prime}:$ $\displaystyle\\{2,0\\},$
$\displaystyle\psi_{2}$ $\displaystyle:=C_{1342}=\psi_{2}^{\prime}:$
$\displaystyle\\{0,0\\},$ $\displaystyle\psi_{3}$
$\displaystyle:=C_{2124}=\psi_{1}^{\prime}:$ $\displaystyle\\{-2,0\\},$
$\displaystyle\psi_{4}$ $\displaystyle:=C_{2424}=\psi_{0}^{\prime}:$
$\displaystyle\\{-4,0\\},$

these are usually known as the NP scalars. Similarly, the 6 components of the
Maxwell tensor, can be encoded in 3 complex scalars:

$\displaystyle\phi_{0}$ $\displaystyle=l^{a}m^{b}F_{ab}:$
$\displaystyle\\{2,0\\}$ (37) $\displaystyle\phi_{1}$
$\displaystyle=\frac{1}{2}\left(l^{a}n^{b}-m^{a}\bar{m}^{b}\right)F_{ab}:$
$\displaystyle\\{0,0\\}$ $\displaystyle\phi_{2}$
$\displaystyle=\bar{m}^{a}n^{b}F_{ab}:$ $\displaystyle\\{-2,0\\}$

with $\phi_{i}^{\prime}=-\phi_{2-i}$.

In certain cases, the algebraic structure of the Weyl tensor greatly
simplifies the formalism. We say that a null vector is tangent to a principal
null direction of the Weyl tensor if it satisfies the following [52]:

$k^{b}k^{c}k_{[e}C_{a]bc[d}k_{f]}=0\,.$ (38)

In general, a spacetime has 4 principal null directions, however, in certain
instances, they may overlap. When they coincide in two pairs of two, we state
that the spacetime is of Petrov type-D. This is the case for Kerr BHs. By
taking $l^{a}$ and $n^{a}$ to be tangent to these null directions, many of the
quantities defined above vanish. In fact, for a vacuum spacetime, we can use
the Goldberg-Sachs theorem to establish that:

$\displaystyle\psi_{0}$ $\displaystyle=\psi_{1}=\psi_{3}=\psi_{4}=0\,,$ (39)
$\displaystyle\kappa$
$\displaystyle=\kappa^{\prime}=\sigma=\sigma^{\prime}=0\,.$

Thus, the only non vanishing GHP scalars are
$\psi_{2},\,\rho,\,\rho^{\prime},\,\tau$ and $\tau^{\prime}$. The Kinnersley
tetrad is an example that satisfies this condition [57, 43]:

$\displaystyle l^{\mu}$
$\displaystyle=\left(\frac{a^{2}+r^{2}}{\Delta_{r}},1,0,\frac{a}{\Delta_{r}}\right)^{\mu},$
(40) $\displaystyle n^{\mu}$
$\displaystyle=\frac{1}{2\Sigma}\left(a^{2}+r^{2},-\Delta,0,a\right)^{\mu},$
$\displaystyle m^{\mu}$
$\displaystyle=\frac{i}{\sqrt{2}\left(r+iax\right)}\left(a\sqrt{\Delta_{x}},0,i\sqrt{\Delta_{x}},\frac{1}{\sqrt{\Delta_{x}}}\right)^{\mu},$
$\displaystyle\bar{m}^{\mu}$
$\displaystyle=\frac{-i}{\sqrt{2}\left(r-iax\right)}\left(a\sqrt{\Delta_{x}},0,-i\sqrt{\Delta_{x}},\frac{1}{\sqrt{\Delta_{x}}}\right)^{\mu}.$

###  II.3 Master wave equation

Teukolsky first proved that the EOM of most fields in a Kerr background can be
described by a single wave equation [3]. Following the approach in [58], we
will express this result in a more covariant way. We start by defining a
generalized d’Alembert operator [58]:

$\hskip 1.42271pt\mathaccent 0{\Box}{\text{\tiny{T}}}\hskip
1.42271pt_{p}=\left(\Theta^{a}+pB^{a}\right)\left(\Theta_{a}+pB_{a}\right)\,,$
(41)

where

$B^{a}=\rho\,n^{a}-\tau\,\bar{m}^{a}\,.$ (42)

and we defined $\Theta_{a}$ in equation (34). When acting on a $\\{0,0\\}$
qunatity, $\Theta_{a}$ reduces to $\nabla_{a}$, thus $\hskip
1.42271pt\mathaccent 0{\Box}{\text{\tiny{T}}}\hskip 1.42271pt_{0}$ is the
standard d’Alembert operator $\Box$. The Teukolsky master equation is a
generalization of this that acts on $\psi^{(s)}$, a spin-s field with weight
$\\{2s,0\\}$ (see table 1). The master Teukolsky equation reads [58]:

$\mathcal{O}_{s}\,\psi^{(s)}:=\frac{1}{2}\left(\hskip 1.42271pt\mathaccent
0{\Box}{\text{\tiny{T}}}\hskip
1.42271pt_{2s}-4s^{2}\psi_{2}\right)\psi^{(s)}=0\,.$ (43)

Using the Kinnersley tetrad (equation (40)), this equation can be written:

$\mathcal{O}_{s}\,\psi^{(s)}=\left((\text{\thorn}-2s\rho-\bar{\rho})(\text{\thorn}^{\prime}-\rho^{\prime})-(\eth-2s\tau-\bar{\tau}^{\prime})(\eth^{\prime}-\tau^{\prime})-(2s^{2}-3s+1)\psi_{2}\right)\psi^{(s)}=0\,.$
(44)

Now, to separate the resulting PDE, we follow the approach in [3]. We make the
Ansatz:

$\psi^{(s)}(t,r,x,\phi)=S(t,\phi)\,\Psi_{s}(r,x)\,,$ (45)

where

$S(t,\phi):=e^{-i\left(\omega t-m\phi\right)}\,,$ (46)

and $m\in\mathbb{Z}$, ensuring periodic boundary conditions are met for
$\phi$. Substituting this in the Teukolsky equation, we obtain:

$\mathcal{O}_{s,\omega}\,\Psi_{s}:=\frac{1}{S}\mathcal{O}_{s}\psi^{(s)}\,=0,$
(47)

where $\mathcal{O}_{s,\omega}$ is an operator that only depends on $(r,x)$.
Finally, taking $\Psi_{s}=R_{s}(r)\Theta_{s}(x)$, we can separate equation
(47) as the sum of a radial and an angular ODE operator [3, 16, 59]:

$\mathcal{O}_{s,\omega}\,\Psi_{s}=\frac{1}{\Sigma(r,x)}\left(L_{r,s}+L_{x,s}\right)\left(\Theta_{s}(x)R_{s}(r)\right)=0\,,$
(48)

where:

$\displaystyle L_{x,s}\Theta_{s}(x)$
$\displaystyle:=\derivative{x}\Delta_{x}\derivative{x}\Theta_{s}(x)+\left(s-\frac{(m+sx)^{2}}{\Delta_{x}}-2a\omega
sx+a^{2}\omega^{2}x^{2}+{}_{s}\mathcal{A}_{\ell m}\right)\Theta_{s}(x)=0\,,$
(49) $\displaystyle L_{r,s}R_{s}(r)$
$\displaystyle:=\Delta_{r}^{-s}\derivative{r}\Delta_{r}^{s+1}\derivative{r}R_{s}(r)+\left(\frac{K^{2}-is\Delta_{r}^{\prime}K}{\Delta_{r}}+4is\omega
r-{}_{s}\Lambda_{lm}\right)R_{s}(r)=0\,,$ (50)

with:

$\displaystyle K(r)$
$\displaystyle=\left(a^{2}+r^{2}\right)\left(\omega-m\frac{a}{a^{2}+r^{2}}\right)\,,$
(51) $\displaystyle{}_{s}\Lambda_{lm}$ $\displaystyle={}_{s}\mathcal{A}_{\ell
m}-a^{2}\omega^{2}+2am\omega\,,$

and ${}_{s}\mathcal{A}_{\ell m}$ is a separation constant. In the
Schwarzschild limit, regular solutions of equation (49) are the spin-weighted
spherical harmonics, with ${}_{s}\mathcal{A}_{\ell m}\rightarrow
l(l+1)-s(s+1)$, see [60]. When we turn on rotation, the equation depends on
$\omega$, and we no longer know of a closed-form solution. Regular solutions
in this regime are known as spin-weighted spheroidal harmonics.

The singular structure of equations (49) and (50) is identical. They both have
two regular singular points at finite values of their domain ($x=\pm 1$ for
(49) and $r=r_{\pm}$ for (50)), and an irregular singular point at infinity.
In fact, through suitable reparametrizations, they can be both expressed as a
confluent Heun equation, see [59, 61].

$s\hskip 3.99994pt$ | $\psi^{(s)}\hskip 3.99994pt$ |
---|---|---
0 | $\Phi\hskip 3.99994pt$ | Scalar field
$\frac{1}{2}$ | $\chi_{0}\hskip 3.99994pt$ | Dirac spinors
$-\frac{1}{2}$ | $\psi_{2}^{-\nicefrac{{1}}{{3}}}\,\chi_{1}\hskip 3.99994pt$
$1$ | $\phi_{0}\hskip 3.99994pt$ | Maxwell fields
$-1$ | $\psi_{2}^{-\nicefrac{{2}}{{3}}}\,\phi_{2}\hskip 3.99994pt$
$2$ | $\psi_{0}{}^{(1)\hskip 2.79996pt}$ | GWs
$-2$ | $\psi_{2}^{-\nicefrac{{4}}{{3}}}\,\psi_{4}{}^{(1)}\hskip 3.99994pt$
Table 1: Physical meaning of the fields that solve equation (43).

###  II.4 Quasinormal Modes

QNMs are the natural basis to encode the late time behaviour of linear fields
in a BH spacetime. Because the system is naturally dissipative, the
corresponding frequencies are complex valued. To formally define these modes,
the standard approach is to introduce hyperboloidal coordinates in a
conformally compactified version of spacetime. QNMs are solutions of a given
wave equation, with defined frequency, that are regular at $\mathcal{H}^{+}$
and $\mathcal{J}^{+}$, see [62, 63, 64]. Physically, this is the statement
that QNMs are waves that are ingoing through $\mathcal{H}^{+}$ and outgoing at
$\mathcal{J}^{+}$. For simplicity, in this paper, we will work with a constant
$t$ slicing of spacetime, and derive the BCs for the radial function that
encode the appropriate ingoing / outgoing behaviour of the full solution. We
will specialize in QNM solutions of the Teukolsky equation, (see (43)).

A Frobenius analysis of the radial EOM near $r=r_{+}$ shows that the radial
function can be decomposed as the sum of an ingoing and an outgoing mode [3,
16]:

$R_{s}(r)=\mathcal{A}(\omega)R_{\text{out},r_{+}}(r)+\mathcal{B}(\omega)R_{\text{in},r_{+}}(r)\,,$
(52)

where:

$\displaystyle R_{\text{out},r_{+}}$
$\displaystyle=\left(r-r_{+}\right)^{i\rho}\left(1+\mathcal{O}(r-r_{+})^{-1}\right)\,,$
(53) $\displaystyle R_{\text{in},r_{+}}$
$\displaystyle=\left(r-r_{+}\right)^{-s-i\rho}\left(1+\mathcal{O}(r-r_{+})^{-1}\right)\,,$

with:

$\rho:=\frac{\omega-m\,\Omega_{+}}{2\,\kappa}\,.$ (54)

Similarly, an asymptotic expansion near $r\rightarrow\infty$ yields a
decomposition into an ingoing and an outgoing mode:

$R_{s}(r)=\mathcal{C}(\omega)R_{\text{out},+\infty}(r)+\mathcal{D}(\omega)R_{\text{in},+\infty}(r)\,,$
(55)

where:

$\displaystyle R_{\text{out},+\infty}(r)$
$\displaystyle=r^{-1-2s+2iM\omega}e^{i\omega
r}\left(1+\mathcal{O}\left(r^{-1}\right)\right)\,,$ (56) $\displaystyle
R_{\text{in},+\infty}(r)$ $\displaystyle=r^{-1-2iM\omega}e^{-i\omega
r}\left(1+\mathcal{O}\left(r^{-1}\right)\right)\,.$

QNMs are defined as solutions of the Teukolsky equation with
$\mathcal{A}(\omega)=\mathcal{D}(\omega)=0$. This defines a boundary value
problem (BVP) with eigenfunction $R_{s}$ and eigenvalue $\omega$.

Numerically, there is a more practical way of setting up the BVP. Define the
compact radial coordinate:

$y:=1-\frac{r_{+}}{r}\,.\\\ $ (57)

At $y=0$ ($r=r_{+}$), and $y=1$ ($r=+\infty$), the BCs corresponding to
ingoing and outgoing modes and the multiplicative factor that transforms one
into the other are singular. Define:

$R_{s}(r(y))=\alpha_{y,s}(y)\widetilde{R}_{s}(y)\,,$ (58)

where $\alpha_{y,s}$ encodes the singular part of the appropriate BCs, [16]:

$\alpha_{y,s}(y(r)):=e^{ir\omega}r^{4iM\omega}\Delta_{r}^{-s-iM\omega}\left(\frac{r-r_{-}}{r-r_{+}}\right)^{\frac{i\left(\frac{M}{r_{+}}\omega-m\Omega\right)}{2\kappa}}\,.$
(59)

Substituting this in equation (50), and requiring smoothness of
$\widetilde{R}_{s}(y)$ in $y\in[0,1]$, the appropriate QNM BCs are
automatically satisfied. As outlined in section A, pseudo-spectral methods can
only encode smooth functions, thus, if we find a pseudo-spectral solution of
the resulting EOM in the interval $y\in[0,1]$, it must describe a QNM.

In modified BL coordinates, the Kerr metric has coordinate singularities at
both $S^{2}$ poles ($x=\pm 1$). This induces apparent singular behaviour on
Teukolsky modes in these regions. As above, we can circumvent this by
explicitly factoring out the problematic component. Take:

$\Theta_{s}(x)=\alpha_{x,s}(x)\widetilde{\Theta}_{s}(x)\,,$ (60)

with:333$\alpha_{x,s}(x)$ is smooth for $x\in[-1,1]$, however, when we combine
this with $e^{-i\,m\,\phi}$ and convert to Cartesian coordinates, we get an
irregular function.

$\alpha_{x,s}(x):=(1-x)^{\frac{\absolutevalue{m+s}}{2}}(1+x)^{\frac{\absolutevalue{m-s}}{2}}\,,$
(61)

QNM solutions will have $\widetilde{\Theta}_{s}(x)$ to be a smooth function
for $x\in[-1,1]$.

This definition of QNMs relies on the separability of the EOM into a radial
and an angular ODE. However, in sections III and IV, we will study EOM that do
not separate. The same set of steps is still valid. First, using a Frobenius
analysis of the PDE near the boundaries of the domain, determine the BCs that
encode the correct ingoing / outgoing behaviour. Then, define a function
$\alpha_{s}(y,x)$ that encodes the corresponding singular behaviour and factor
it out from the QNM solution. Finally, solve the resulting EOM requiring
smoothness of the solution for $[y,x]\in[0,1]\crossproduct[-1,1]$. Concretely,
for Teukolsky modes, the joint singular behaviour is:

$\alpha_{s}(y,x)=\alpha_{y,s}(y)\alpha_{x,s}(x)\,.$ (62)

Define:

$\widetilde{\Psi}_{s}:=\frac{1}{\alpha_{s}}\Psi_{s}\,.$ (63)

Substituting this in equation (47), we get:

$0=\widetilde{\mathcal{O}}_{s,\omega}\,\widetilde{\Psi}_{s}(y,x):=\frac{1}{\alpha_{s}}\mathcal{O}_{s,\omega}\,\alpha_{s}\widetilde{\Psi}_{s}(y,x)\,.$
(64)

QNMs are solution of the full PDE that are smooth for
$[y,x]\in[0,1]\crossproduct[-1,1]$.

###  II.5 Parity transformations

Denote $\mathbb{P}$ the parity transformation around $r=0$. In modified BL
coordinates, $\mathbb{P}$ takes $x\rightarrow-x$ and
$\phi\rightarrow\pi+\phi$. The action (6) is invariant under parity, and so is
the Kerr metric. For the Kinnersley tetrad,see (40), we can prove:

$\displaystyle\mathbb{P}l^{a}\partial_{a}$ $\displaystyle=l^{a}\partial_{a},$
(65) $\displaystyle\mathbb{P}n^{a}\partial_{a}$
$\displaystyle=n^{a}\partial_{a},$ $\displaystyle\mathbb{P}m^{a}\partial_{a}$
$\displaystyle=-\bar{m}^{a}\partial_{a},$
$\displaystyle\mathbb{P}\bar{m}^{a}\partial_{a}$
$\displaystyle=-m^{a}\partial_{a}.$

The vector potential $A_{a}$, corresponding field strength $F_{ab}$, and the
spacetime metric are all real quantities. However, because of the complex
nature of $m$ and $\bar{m}$, when we project these tensors to the null tetrad,
the resulting scalars are complex-valued. Notwithstanding, this is merely an
artifact of the representation. In particular, the complex conjugation
operator simply maps $m\rightarrow\bar{m}$ and $\bar{m}\rightarrow m$. This
motivates the definition of the conjugate-parity-transform operator
$\mathcal{P}$. For some scalar $\eta$, we have:

$\mathcal{P}\eta:=\overline{\mathbb{P}\eta}\,.$ (66)

As an example, consider the action in the Teukolsky field $\psi^{(}s)$. Using
equation (45), we have that:

$\displaystyle\mathcal{P}\psi^{(s)}(t,r,x,\phi)$
$\displaystyle=\mathcal{P}S(t,\phi)\,\Psi_{s}(r,x)$ (67)
$\displaystyle=(-1)^{m}\bar{S}(t,\phi)\bar{\Psi}_{s}(r,-x)\,,$

where the $(-1)^{m}$ contribution comes from taking $\phi\rightarrow\pi+\phi$.

The transformation rule in equation (65) is tensorial. Thus, if $\eta$ is
constructed by taking covariant derivatives and tensor products of tetrad
quantities, we have:

$\mathcal{P}\eta=(-1)^{s}\eta\,,$ (68)

where $s$ is the number of $m^{a}$ and $\bar{m}^{a}$ instances in the
definition of $\eta$. As an example, take $\rho$ and $\tau$. Using equation
(25), we have:

$\mathcal{P}\rho=(-1)^{2}\rho\qquad\mathcal{P}\tau=(-1)^{1}\tau\,.$ (69)

By definition, $\eta$ is a GHP covariant quantity. Thus, $s$ is simply its
spin weight: $s=(p-q)/2$.

Now, consider $\mathcal{G}$, a differential operator constructed with tensor
products and covariant derivatives of tetrad quantities, acting on some scalar
$f$. Following the same reasoning, we have:

$\mathcal{P}\,\left(\mathcal{G}\,f\right)=(-1)^{s}\mathcal{G}\,\mathcal{P}\,f\,.$
(70)

where $s$ is the spin-weight of $\mathcal{G}$. Because Kerr is invariant under
parity, this rule also extends to the case where $\mathcal{G}$ contains pure
gravitational scalars, e.g. $\psi_{2}$.

In particular, the Teukolsky operator $\mathcal{O}_{s}$ has spin-weight $0$,
thus:

$\left[\mathcal{P},\mathcal{O}_{s}\right]=0,$ (71)

where $[\cdot,\cdot]$ denotes the commutator operator. Standard results tell
us that if $\psi^{(s)}$ is in the kernel of $\mathcal{O}_{s}$, then so is
$\mathcal{P}\psi^{(s)}$. Using equation (45), we have that:

$\displaystyle\mathcal{O}_{s}\,S\,\Psi_{s}=0\implies$ (72)
$\displaystyle\mathcal{O}_{s}\,(-1)^{m}\bar{S}\,\mathbb{P}\bar{\Psi}_{s}=0\,.$

This proves that if $\left(m,\omega\right)$ is in the QNM spectrum of
$\mathcal{O}_{s}$, then so is $\left(-m,-\bar{\omega}\right)$. The parity
invariance of the Kerr metric leads to a degeneracy between these two families
of QNMs.

The action of $\mathbb{P}$ on the Kinnersley tetrad, see (65), does not
generalize to all the equivalent tetrads. In general, if we generate a new
null tetrad with a GHP transformation, see (21), the identities in (65) are no
longer respected. To circumvent this, we can restrict the allowed GHP
transformations, such that the resulting tetrads obey (65). Under a generic
GHP transformation, we have:

$\displaystyle\mathbb{P}l^{a}\partial_{a}$
$\displaystyle\rightarrow\mathbb{P}\left(\lambda\bar{\lambda}\right)\mathbb{P}l^{a}\partial_{a},$
(73) $\displaystyle\mathbb{P}n^{a}\partial_{a}$
$\displaystyle\rightarrow\mathbb{P}\left(\frac{1}{\lambda\bar{\lambda}}\right)\mathbb{P}n^{a}\partial_{a},$
$\displaystyle\mathbb{P}m^{a}\partial_{a}$
$\displaystyle\rightarrow\mathbb{P}\left(\frac{\lambda}{\bar{\lambda}}\right)\mathbb{P}m^{a}\partial_{a},$
$\displaystyle\mathbb{P}\bar{m}^{a}\partial_{a}$
$\displaystyle\rightarrow\mathbb{P}\left(\frac{\bar{\lambda}}{\lambda}\right)\mathbb{P}\bar{m}^{a}\partial_{a}.$

Thus, to preserve (65), we must have:

$\displaystyle\mathbb{P}\lambda\bar{\lambda}$
$\displaystyle=\lambda\bar{\lambda}\,,$ (74)
$\displaystyle\mathbb{P}\left(\frac{\lambda}{\bar{\lambda}}\right)$
$\displaystyle=\frac{\bar{\lambda}}{\lambda}\,.$

The absolute value and phase of $\lambda$ must be parity even and odd
respectively. Take $\eta$ to be a GHP scalar with boost weight $r$ and spin
weight $s$. Restricting to GHP transformations that respect (74), we get:

$\displaystyle\mathbb{P}\eta$
$\displaystyle\rightarrow\left(\lambda\bar{\lambda}\right)^{r}\left(\frac{\lambda}{\bar{\lambda}}\right)^{-s}\mathbb{P}\eta\,,$
(75) $\displaystyle\mathcal{P}\eta$
$\displaystyle\rightarrow\left(\lambda\bar{\lambda}\right)^{r}\left(\frac{\lambda}{\bar{\lambda}}\right)^{s}\mathcal{P}\eta\,.$

The action $\mathbb{P}$ is akin to complex conjugation. It takes a $\\{p,q\\}$
GHP scalar to a $\\{q,p\\}$ scalar. Consequently, $\mathcal{P}$ preserves the
GHP weight. This fact is particularly useful, as it allows us to take linear
combinations of $\eta$ and $\mathcal{P}\eta$ without breaking GHP covariance.

###  II.6 Hodge duals and p forms

In this section, we will review some key concepts regarding $p$ forms, and
their hodge duals. The hodge dual operator, $\star$, is a linear map from the
space $p$ forms to the space of $d-p$ forms. We will follow the convention:

$\left(\star\omega\right)_{a_{1}\cdots
a_{d-p}}=\frac{1}{p!}\varepsilon_{a_{1}\cdots a_{d-p}b_{1}\cdots
b_{p}}\omega^{b_{1}\cdots b_{p}}\,,$ (76)

where $\varepsilon$ is the Levi-Civita tensor. We have
$\varepsilon_{a_{1}\cdots a_{d}}=\varepsilon_{[a_{1}\cdots a_{d}]}$ and
$\varepsilon_{123\cdots d}=\sqrt{-g}$. For a Lorentzian metric, we have:

$\star\star\omega=-(-1)^{p\left(d-p\right)}\omega\,.$ (77)

Finally, we can define the differential and co-differential operators:

$\displaystyle\left(\differential\omega\right)_{a_{1}\cdots a_{p+1}}$
$\displaystyle:=(p+1)\nabla_{[a_{1}}\omega_{a_{2}\cdots a_{p+1}]}\,,$ (78)
$\displaystyle\left(\differential^{\dagger}\omega\right)_{a_{1}\cdots
a_{p-1}}$
$\displaystyle:=\left(\star\differential\star\omega\right)_{a_{1}\cdots
a_{p-1}}$
$\displaystyle=-(-1)^{p\left(d-p\right)}\nabla^{b}\omega_{a_{1}\cdots
a_{p-1}b}\,.$

The Weyl tensor is not fully antisymmetric, and thus cannot be regarded as a
$4$ form. However, it is antisymmetric in the first two and the last two
indices, and so, we can take their dual. In 4-dimensions, we define the left
and right dual of $\mathcal{C}$:

$\displaystyle\left(\star\mathcal{C}\right)_{abcd}=\frac{1}{2!}\varepsilon_{abef}\mathcal{C}^{ef}{}_{cd}\,,$
(79)
$\displaystyle\left(\mathcal{C}\star\right)_{abcd}=\frac{1}{2!}\mathcal{C}_{ab}{}^{ef}\varepsilon_{efcd}\,.$

Now, taking the left and right dual at the same time, we get:

$\displaystyle\left(\star\mathcal{C}\star\right)^{ab}{}_{cd}$
$\displaystyle=\frac{1}{4}\varepsilon^{abef}\varepsilon_{cdgh}\mathcal{C}_{ef}{}^{gh}$
(80)
$\displaystyle=-\frac{4!}{4}\delta^{a}_{[c}\delta^{b}_{d}\delta^{e}_{g}\delta^{f}_{h]}\,\mathcal{C}_{ef}{}^{gh}$
$\displaystyle=-\frac{4!}{4\cdot
3!}\delta^{e}_{[c}\delta^{f}_{d]}\delta^{a}_{[g}\delta^{b}_{h]}\,\mathcal{C}_{ef}{}^{gh}$
$\displaystyle=-\mathcal{C}^{ab}{}_{cd}\,.$

In the second step we used the standard result for the contraction of two
instances of $\varepsilon$, and in the the third step, we used the fact that
$\mathcal{C}$ is traceless, i.e. $\mathcal{C}^{c}{}_{acb}=0$. Using
$\star^{2}=-1$ (for a 2 form in 4 dimensions), we get:

$\star\mathcal{C}=\mathcal{C}\star\,,$ (81)

The left and right duals are equivalent.

There are multiple ways to fix the degeneracy between the two electromagnetic
/ gravitational degrees of freedom. Usually, it is helpful to decompose the
relevant quantities in a basis of eigenstates of a given operator. If we pick
that operator to be $\star$, we decompose the electromagnetic and Weyl tensor
into their self-dual and anti-self-dual parts:

$\displaystyle F^{\pm}$ $\displaystyle=\frac{1}{2}\left(F\mp i\star
F\right)\,,$ (82) $\displaystyle\mathcal{C}^{\pm}$
$\displaystyle=\frac{1}{2}\left(\mathcal{C}\mp i\star\mathcal{C}\right)\,,$

As intended, $F^{\pm}$ and $\mathcal{C}^{\pm}$ are $\star$ eigenstates:

$\star F^{\pm}=\pm iF^{\pm}\,,\qquad\star\mathcal{C}^{\pm}=\pm
i\mathcal{C}^{\pm}\,.$ (83)

When projected to the Kinnersley tetrad (40), the Levi-Civita tensor takes a
simple form:

$\displaystyle\varepsilon_{(\alpha)(\beta)(\gamma)(\delta)}$
$\displaystyle=\varepsilon_{abcd}e^{a}_{(\alpha)}e^{b}_{(\beta)}e^{c}_{(\gamma)}e^{d}_{(\delta)}$
(84) $\displaystyle=-i\,E_{(\alpha)(\beta)(\gamma)(\delta)}$

where $E$ is fully antisymmetric and $E_{1234}$. Thus, we get
$\varepsilon_{(1)(2)(3)(4)}=-i$. Dualized quantities, carry an extra factor of
$i$, and thus, complex conjugation will add an extra $-$ contribution.
Concretely, consider the dualized GB term:

$\displaystyle\widetilde{\mathcal{G}}$
$\displaystyle=\left(\star\mathcal{C}\right)^{abcd}\mathcal{C}_{abcd}$ (85)
$\displaystyle=24i\left(\bar{\psi}_{2}^{2}-\psi_{2}^{2}\right)\,.$

$\widetilde{\mathcal{G}}$ can be understood as a product of $i$ and a spin
weight $0$ GHP scalar that obeys the transformation rules detailed in section
II.5. Thus, under the conjugate-parity-transform, we have:

$\mathcal{P}\widetilde{\mathcal{G}}=-\widetilde{\mathcal{G}}\,.$ (86)

As a rule of thumb, quantities that are dualized an odd amount of times are
parity odd.

###  II.7 The vacuum Maxwell equations

In the GHP formalism, we represent the electromagnetic field with 3 complex
scalar fields. However, we know that electromagnetism is fully specified by
two independent polarizations. Thus, two of the GHP scalars must be redundant.
In [65], Wald derived an approach to fully reconstruct the electromagnetic
tensor $F_{ab}$ from a single complex scalar. This is commonly known in the
literature as the Hertz potential approach. In this section, we will start by
deriving the Teukolsky equation for $\phi_{0}$, and then explain how to
generate $F_{ab}$ from it.

The vacuum Maxwell equations take the form:

$\displaystyle\nabla^{a}F_{ab}=0,$ (87) $\displaystyle\differential F=0.$

Here, the second equation, commonly known as the Bianchi identity, implies
that $F_{ab}$ must be a closed $2$ form. Because $F$ is a real field, we can
use the results in section II.6, to combine the two equations into a single
complex equation:

$M_{a}:=i\left(\star\differential F^{-}\right)_{a}=0\,,$ (88)

where $F^{-}$ was defined in (82). The, first and second equations in (87),
are encoded by the real and imaginary parts of (88) respectively. In the
tetrad basis, using the GHP formalism, $M$ takes a simple form [53]:

$M_{(\alpha)}=\begin{pmatrix}\left(\eth^{\prime}-\tau^{\prime}\right)\phi_{0}-\left(\text{\thorn}-2\rho\right)\phi_{1}\\\
\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)\phi_{1}-\left(\eth-\tau\right)\phi_{2}\\\
\left(\text{\thorn}^{\prime}-\rho^{\prime}\right)\phi_{0}-\left(\eth-2\tau\right)\phi_{1}\\\
\left(\eth^{\prime}-2\tau^{\prime}\right)\phi_{1}-\left(\text{\thorn}-\rho\right)\phi_{2}\end{pmatrix}_{(\alpha)}\,.$
(89)

The corresponding GHP weights are:

$\\{1,1\\},\quad\\{-1,-1\\},\quad\\{1,-1\\},\quad\\{-1,1\\}\,.$ (90)

If we increase the spin-weight of $M_{l}$ and the boost-weight of $M_{m}$ by
1, the two components will have the same GHP weights and can be combined. That
is precisely the role of the operator $\zeta$ introduced by Teukolsky in [3]:

$\zeta^{a}:=\left(\eth-2\tau+\bar{\tau}^{\prime}\right)l^{a}-\left(\text{\thorn}-2\rho+\bar{\rho}\right)m^{a}\,.$
(91)

Using the GHP relations in [53], we can obtain the Teukolsky equation for
$s=1$ (see equation (44)):

$\zeta^{a}M_{a}=\mathcal{O}_{1}\phi_{0}=0\,.$ (92)

$F$ is a closed form, thus, locally, there is a form $A$, such that
$F=\differential A$. This is known as the vector potential. We can express the
components of $F$ directly in terms of the $A$. Using the definition of
$\phi_{i}$ (equation (37)), and treating $A^{a}$ as a $\\{0,0\\}$ GHP
covariant vector field, we have:

$\displaystyle\phi_{0}$
$\displaystyle=\tau_{0}{}^{a}A_{a}:=\left(m^{a}\text{\thorn}-l^{a}\eth\right)A_{a}\,,$
(93) $\displaystyle\phi_{1}$
$\displaystyle=\tau_{1}{}^{a}A_{a}:=\frac{1}{2}\left(n^{a}\text{\thorn}-l^{a}\text{\thorn}^{\prime}+m^{a}\eth^{\prime}-\bar{m}^{a}\eth\right)A_{a}\,,$
$\displaystyle\phi_{2}$
$\displaystyle=\tau_{2}{}^{a}A_{a}:=\left(n^{a}\eth^{\prime}-\bar{m}^{a}\text{\thorn}^{\prime}\right)A_{a}\,.$

Thus, in terms of $A$, the Teukolsky equation is:

$\mathcal{O}_{1}\,\tau_{0}\,A=0\,.$ (94)

On the other hand, working directly in an abstract basis, we can obtain
$M_{a}$ in terms of $A_{a}$. We have:

$M_{a}=\left(\mathcal{E}A\right)_{a}:=\nabla_{b}\nabla^{b}A_{a}-\nabla_{b}\nabla_{a}A^{b}\\\
$ (95)

Acting with $\zeta^{a}$ on the left, we must recover the Teukolsky equation:

$\mathcal{O}_{1}\phi_{0}=\zeta\,\mathcal{E}\,A\,.$ (96)

Thus, we have the operator identity:

$\mathcal{O}_{1}\tau_{0}=\zeta\mathcal{E}\,.$ (97)

The key insight comes from taking the transpose444Here the transpose refers to
the classical transpose of a differential operator, no complex conjugation is
needed. See [65]for details. of this identity. We have that:

$\tau_{0}{}^{T}\mathcal{O}_{1}{}^{T}=\mathcal{E}^{T}\zeta^{T}\,.$ (98)

Now, we can prove that $\mathcal{E}^{T}=\mathcal{E}$, thus, for a mode
$\psi_{H}$ such that:

$\mathcal{O}_{1}{}^{T}\psi_{H}=0\,,$ (99)

we have:

$\mathcal{E}\left(\zeta^{T}\psi_{H}\right)=0\,.$ (100)

We proved that $\hat{A}_{a}=\zeta_{a}{}^{T}\psi_{H}$, solves the vacuum
Maxwell equations. Thus, by solving equation (99), we can generate all the
components of $F_{ab}$. This is not yet the complete story. In general,
$\psi_{H}$ is a complex scalar, thus, $\hat{A}_{a}$ is a complex vector.
Nevertheless, the reality of $A_{a}$ and consequently $F_{ab}$ is a crucial
assumption in the derivation of equation (88). $\mathcal{E}$ is a real linear
operator, thus if $\hat{A}$ is in its kernel, then so must be
$\overline{\hat{A}}$. Hence, for any mode $\psi_{H}$ that solves (99), we can
construct a real solution to Maxwell’s equation given by:

$A_{a}:=\mathcal{R}e(\zeta_{a}{}^{T}\psi_{H}).$ (101)

In the literature, $\psi_{H}$ is usually known as the Hertz potential.

From equation (43), we can deduce that
$\mathcal{O}_{s}{}^{T}=\mathcal{O}_{-s}$. Thus, $\psi_{H}$ is a solution to
the master Teukolsky equation with $s=-1$. In fact, this is a re-
parametrization of the EOM for $\phi_{2}$ [58]:

$\mathcal{O}_{-1}\left(\psi_{2}^{-2/3}\phi_{2}\right)=0$ (102)

Thus, the Hertz potential is proportional to the $\phi_{2}$ component of a
solution to the Maxwell equations, and we can deduce $\psi_{H}$ has GHP weight
$\\{-2,0\\}$. We can prove that any QNM of $\phi_{2}$ will generate a QNM for
$A^{a}$ with the same frequency, and consequently, all QNMs of $\phi_{2}$ must
be QNMs of $\phi_{0}$. The Hertz potential derivation can be reversed, and we
can generate a vector potential $A^{a}$ from any $\phi_{0}$. Thus, all
$\phi_{0}$ QNMs must also be QNMs of $\phi_{2}$. Hence, $\phi_{0}$ and
$\phi_{2}$ must be isospectral, as seen in [43].

Notwithstanding, the $\phi_{2}$ proportional to the Hertz potential, is not
the same as the component of the electromagnetic tensor generated by (101).
Instead, we have to use the operators $\tau_{i}$ defined in (93). Applying the
GHP identities in [53] we get [65]:

$\displaystyle\phi_{0}$
$\displaystyle=\tau_{0}{}^{a}A_{a}=\frac{1}{2}\text{\thorn}\text{\thorn}\bar{\psi}_{H}$
(103) $\displaystyle\phi_{1}$
$\displaystyle=\tau_{1}{}^{a}A_{a}=\frac{1}{2}\left(\text{\thorn}\,\eth^{\prime}+\tau^{\prime}\,\text{\thorn}\right)\bar{\psi}_{H}$
$\displaystyle\phi_{2}$
$\displaystyle=\tau_{2}{}^{a}A_{a}=\frac{1}{2}\eth^{\prime}\eth^{\prime}\bar{\psi}_{H}$

Because we required $A$ to be a real vector, $\bar{\phi}_{i}$ can be obtained
by simply taking the complex conjugate of these expressions. If we worked with
complex $A$, $\phi_{i}$ and $\bar{\phi}_{i}$ would be independent variables.
Conveniently, $\phi_{i}$ depends exclusively on $\bar{\psi}_{H}$. 555This
coincidence is merely a consequence of the definition of $\phi_{i}$. If the
projection of $F$ into the tetrad basis was done differently, we would likely
have a mixing of $\psi_{H}$ and $\bar{\psi}_{H}$. Thus, we can encode equation
(103), in operators $\mathcal{T}_{i}$ such that:

$\mathcal{T}_{i}\bar{\psi}_{H}:=\phi_{i}$ (104)

##  III Corrections to scalar modes

In this section, we will compute the EFT corrections to QNM frequencies of a
massless scalar, when considering a Kerr background. We will mostly focus in
the parity even corrections, as seen in the action (1), however, in section
III.2, we will briefly outline how to adapt our methods to consider parity odd
corrections.

###  III.1 Parity even equations of motion

The EOM considering parity even EFT corrections, were obtained in (4). For a
Kerr background, using the results in section II.3, they reduce to:

$\left(\mathcal{O}_{0}-L^{2}\lambda\,\mathcal{G}\right)\Phi=0\,,$ (105)

where $\mathcal{O}_{0}$ is the $s=0$ Teukolsky operator and $\mathcal{G}$ is
the GB term, defined in (2). Because Kerr is a vacuum, type-D spacetime, we
get:

$\mathcal{G}=24\left(\psi_{2}^{2}+\bar{\psi}_{2}^{2}\right)\,.$ (106)

Following the discussion in section II.5, it is clear that the conjugate-
parity-transform operator, $\mathcal{P}$, commutes with all terms in (105):

$[\mathcal{P},\mathcal{O}_{0}]=[\mathcal{P},\mathcal{G}]=0\,,$ (107)

thus, if $\Phi$ is a solution to the EOM, then so is $\mathcal{P}\Phi$. As
argued in the end of section II.5, Assuming $\Phi$ is proportional to
$S(t,\phi)=\exp(-i\omega t+im\phi)$, we get that if $(m,\omega)$ define a QNM,
then so does $(-m,-\bar{\omega})$. This relates QNMs whose frequency have
positive real part with the ons with negative real part, and thus, we can
restrict to QNM frequencies with positive real part.

In the Kinnersley tetrad basis,

$\psi_{2}=-\frac{M}{\left(r-i\,ax\right)^{3}}$ (108)

Thus, the EOM are, in general, not separable. Notwithstanding, we can still
follow the procedure in section II.4 to define the QNMs. First, as seen in
(45), we factor out the $(t,\phi)$ dependence, by making the Ansatz:

$\Phi(t,r,x,\phi)=S(t,\phi)\,\Psi(r,x)$ (109)

We get:

$\mathcal{O}_{0,\omega}\,\Psi(r,x)=\lambda L^{2}\,\mathcal{G}\Psi(r,x)$ (110)

where $\mathcal{O}_{0,\omega}$ was defined in equation (47). Now, we must
enforce BCs corresponding to outgoing behavior at $\mathcal{J}^{+}$ and
ingoing at $\mathcal{H}^{+}$. To do so, we will replace $r$ with the
compactified radial coordinate $y=1-r/r_{+}$. In the case $\lambda=0$, as
argued in section II.4, we can set the BCs by factoring out $\alpha_{0}(r,x)$
(see equation (62)) and requiring the remainder to be smooth for
$(y,x)\in[-1,1]\crossproduct[0,1]$. $\mathcal{G}$ is regular at $r_{+}$ and at
$x=\pm 1$. Furthermore, it is subleading with respect to
$\mathcal{O}_{0,\omega}$ at $r\rightarrow\infty$. This implies that the
leading order of a Frobenius expansion at the domain boundaries must be
independent of $\lambda$. Thus, QNMs of the corrected equation obey the same
BCs as the QNMs in the uncorrected case. Define:

$\widetilde{\Psi}(y,x):=\frac{1}{\alpha_{0}(y,x)}\Psi(r(y),x)\,.$ (111)

Replacing this in equation (110), we get:

$\widetilde{\mathcal{O}}_{0,\omega}\,\widetilde{\Psi}(y,x)=L^{2}\lambda\,\mathcal{G}\,\widetilde{\Psi}(y,x)\,.$
(112)

QNMs are solutions of this equation, that are smooth for
$(y,x)\in[-1,1]\crossproduct[0,1]$.

Parameters | $a/M$
---|---
$\ell$ | $m$ | Component | $0$ | $0.25$ | $0.5$ | $0.75$ | $0.95$ | $0.99$ | $0.999$ | $0.9999$ | $0.999975$
$0$ | $0$ | $M\mathcal{R}e(\omega_{0})$ | $0.110455$ | $0.110943$ | $0.112381$ | $0.114322$ | $0.111992$ | $0.110444$ | $0.110265$ | $0.110247$ | $0.110246$
$M\mathcal{I}m(\omega_{0})$ | $-0.104896$ | $-0.104289$ | $-0.102183$ | $-0.097298$ | $-0.089559$ | $-0.089493$ | $-0.089439$ | $-0.089434$ | $-0.089433$
$100\,M\mathcal{R}e(\omega_{1})$ | $14.056932$ | $13.842350$ | $13.018046$ | $10.639102$ | $5.154504$ | $6.752369$ | $7.113310$ | $7.148748$ | $7.151690$
$100\,M\mathcal{I}m(\omega_{1})$ | $-2.446792$ | $-2.302637$ | $-1.884293$ | $-1.484410$ | $-5.668677$ | $-8.548625$ | $-8.774436$ | $-8.796021$ | $-8.797808$
$1$ | $-1$ | $M\mathcal{R}e(\omega_{0})$ | $0.292936$ | $0.275671$ | $0.261572$ | $0.249698$ | $0.241372$ | $0.239810$ | $0.239462$ | $0.239428$ | $0.239425$
$M\mathcal{I}m(\omega_{0})$ | $-0.097660$ | $-0.097362$ | $-0.096505$ | $-0.095278$ | $-0.094124$ | $-0.093882$ | $-0.093828$ | $-0.093822$ | $-0.093822$
$100\,M\mathcal{R}e(\omega_{1})$ | $4.631840$ | $3.835932$ | $3.130117$ | $2.487689$ | $2.055239$ | $1.982260$ | $1.966471$ | $1.964905$ | $1.964775$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.066211$ | $0.362974$ | $0.571957$ | $0.660057$ | $0.625281$ | $0.610799$ | $0.607378$ | $0.607033$ | $0.607005$
$1$ | $1$ | $M\mathcal{R}e(\omega_{0})$ | $0.292936$ | $0.314934$ | $0.344753$ | $0.390378$ | $0.462261$ | $0.493423$ | $0.503344$ | $0.500896$ | $0.500198$
$M\mathcal{I}m(\omega_{0})$ | $-0.097660$ | $-0.097005$ | $-0.094395$ | $-0.086394$ | $-0.060091$ | $-0.036712$ | $-0.017519$ | $-0.007904$ | $-0.004088$
$100\,M\mathcal{R}e(\omega_{1})$ | $4.631840$ | $5.590809$ | $6.833001$ | $8.512134$ | $9.133060$ | $6.466215$ | $2.063483$ | $-0.540870$ | $-0.129697$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.066211$ | $-0.312890$ | $-0.793476$ | $-1.424866$ | $-2.149343$ | $-2.470718$ | $-2.627015$ | $-0.405923$ | $-0.009214$
$2$ | $-2$ | $M\mathcal{R}e(\omega_{0})$ | $0.483644$ | $0.450060$ | $0.422751$ | $0.399883$ | $0.384002$ | $0.381041$ | $0.380384$ | $0.380318$ | $0.380313$
$M\mathcal{I}m(\omega_{0})$ | $-0.096759$ | $-0.096432$ | $-0.095616$ | $-0.094525$ | $-0.093536$ | $-0.093330$ | $-0.093283$ | $-0.093279$ | $-0.093278$
$100\,M\mathcal{R}e(\omega_{1})$ | $2.545392$ | $1.945246$ | $1.494968$ | $1.151127$ | $0.937867$ | $0.900983$ | $0.892928$ | $0.892128$ | $0.892061$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.160232$ | $0.266663$ | $0.315865$ | $0.321760$ | $0.304842$ | $0.300007$ | $0.298872$ | $0.298758$ | $0.298748$
$2$ | $2$ | $M\mathcal{R}e(\omega_{0})$ | $0.483644$ | $0.526740$ | $0.585990$ | $0.679504$ | $0.840982$ | $0.928028$ | $0.978023$ | $0.993234$ | $0.996634$
$M\mathcal{I}m(\omega_{0})$ | $-0.096759$ | $-0.096148$ | $-0.093494$ | $-0.085026$ | $-0.056471$ | $-0.031063$ | $-0.010980$ | $-0.003525$ | $-0.001761$
$100\,M\mathcal{R}e(\omega_{1})$ | $2.545392$ | $3.384658$ | $4.653470$ | $6.805746$ | $9.331381$ | $7.636031$ | $3.618486$ | $1.324989$ | $0.675931$
$100\,M\mathcal{I}m(\omega_{1})$ | $0.160232$ | $-0.026160$ | $-0.338287$ | $-0.872204$ | $-1.415846$ | $-1.032269$ | $-0.322180$ | $-0.038936$ | $-0.000052$
Table 2: Numerical values for the QNM frequencies and corresponding EFT
correction of scalar modes. All the included digits are significant.

(a)

(b)

Figure 4: In this plot we show that our numerical approach converges
exponentially with increasing grid size. Here, $N_{x}$, $N_{r}$ denote the
sizes of the angular and radial grids respectively. We computed the EFT
corrections to scalar and electromagnetic QNMs of BHs with $a/M=0.999975$, and
$\ell=m=1$, for multiple values of $N_{x}$ and $N_{r}$. Then, we compared them
with a reference $\omega_{1}$ obtained with the largest grid, plotting the
difference. On the left panel, we fixed $N_{r}=250$ and varied $N_{x}$,
whereas on the right panel we fixed $N_{x}=75$, and varied $N_{r}$. The
electromagnetic case was particularly difficult to compute, to achieve
accurate results we had to obtain $\omega_{0}$, with $70$ digits of precision.

(a)

(b)

Figure 5: Here, we show a parametric plot of the real and imaginary parts of
the scalar QNM frequencies as $a$ is varied. On the left, we have the
background QNM frequencies, and on the right the EFT correction. To produce
this image, we computed the QNM frequencies of BHs with $a/M\in(0,0.999975)$,
by solving equation (116b). The plot markers represent the QNM frequencies
obtained for the Schwarzschild limit using an independent non-perturbative
method. We find perfect agreement between the two results. $\ell=m\neq 0$
modes vary very quickly near extremality, with $M\omega\rightarrow\ell/2$ as
$a\rightarrow M$. If we truncated the spin at $a/M\sim 0.99$, we would lose
roughly the last third of the parametric curves.

(a)

(b)

Figure 6: In this image we explore the scalar QNM frequencies of $\ell=m\neq
0$ modes of BHs very close to extremality. On the left panel, we show the
background frequency, whereas the right has the EFT correction. We establish
that, as $a\rightarrow M$, $M\omega_{0}\rightarrow\ell/2$, agreeing with the
prediction in [38]. Similarly, in this limit, $M\omega_{1}$ converges to $0$.
Notice that the near-extremal behaviour of these modes has some structure. In
particular, $M\mathcal{R}e(\omega_{0}^{\ell=1})$ crosses $\ell/2$ at
$a/M\approx 0.992$, $M\mathcal{R}e(\omega_{1}^{\ell=1})$ crosses $0$ at
$a/M\approx 0.9993$, and $M\mathcal{I}m(\omega_{1}^{\ell=2})$ crosses $0$ at
$a/M\approx 0.99994$.  Figure 7: Here, we show $M\mathcal{I}m(\omega_{0})$ for
scalar modes in the Schwarzschild limit as a function of $\ell$. The dark line
represents the Eikonal limit prediction. Contrary to gravitational QNMs,
$\mathcal{I}m(\omega_{1})$ increases with increasing $\ell$, approaching the
Eikonal limit frequency from below. This indicates that the longest-lived
modes will have $\ell\rightarrow\infty$.

In the Schwarzschild limit, the RHS reduces to
$\left(48L^{2}M^{2}\lambda/(r_{+}(1-y))^{6}\right)\widetilde{\Psi}(y,x)$, thus
the EOM separate into an angular and a radial ODEs. This is expected due to
the spherical symmetry of the problem. The angular part can be solved
analytically, yielding the familiar spherical harmonics $Y_{\ell m}$. The
radial part is an eigenvalue problem for $\omega$, with $(L/M)^{2}\lambda$ and
$\ell$ as free parameters. This can be solved directly with the pseudospectral
methods discussed in the appendix A. The discussion is more involved when we
consider the spinning case ($a\neq 0$). While the left hand side (LHS) still
separates into an angular and a radial operator, this is no longer the case
for the RHS. To make progress, we will use the validity condition for the EFT
approximation. The action (1) is at most quadratic in $L/M$, thus our results
are only valid to that order. Assuming $\lambda=\mathcal{O}(1)$, the relevant
perturbation parameter is:

$\widetilde{\lambda}:=L^{2}/M^{2}\lambda\,.$ (113)

Our results are only valid to leading order in $\widetilde{\lambda}$.

We can solve the EOM using standard quantum mechanical perturbation theory. We
make the following Ansatz:

$\displaystyle\widetilde{\Psi}(y,x)$
$\displaystyle=\widetilde{\Psi}_{(0)}(y,x)+\widetilde{\lambda}\,\widetilde{\Psi}_{(1)}(y,x)+\mathcal{O}(\widetilde{\lambda}^{2})\,,$
(114) $\displaystyle\omega$
$\displaystyle=\omega_{0}+\widetilde{\lambda}\,\omega_{1}+\mathcal{O}(\widetilde{\lambda}^{2})\,.$
(115)

Grouping together same powers of $\widetilde{\lambda}$, and neglecting higher
order terms the EOM reduce to:

$\displaystyle\widetilde{\mathcal{O}}_{0,\omega_{0}}\,\widetilde{\Psi}_{(0)}=0\,,$
(116a)
$\displaystyle\widetilde{\mathcal{O}}_{0,\omega_{0}}\,\widetilde{\Psi}_{(1)}=-\omega_{1}\partial_{\omega_{0}}\widetilde{\mathcal{O}}_{0,\omega_{0}}\widetilde{\Psi}_{(0)}+M^{2}\mathcal{G}\,\widetilde{\Psi}_{(0)}\,.$
(116b)

Equation (116a) is simply the $s=0$ of (64). Consequently, its solutions will
be the QNMs and frequencies of a massless scalar field. Henceforth, we will
denote these, the background modes. Substituting these in equation (116b), the
EOM are fully specified, and we get an eigenvalue problem for
$(\widetilde{\Psi}_{(1)},\omega_{1})$.

Apart from $(\widetilde{\Psi}_{(1)},\omega_{1})$, equation (116b) is fully
specified by the background quantities. Thus, after solving the background
equation, the EOM for the correction reduce to an eigenvalue problem.

###  III.2 On parity odd EFT corrections

To achieve full generality, we must consider parity odd EFT corrections. For a
vacuum spacetime, under the assumptions we made when deriving the parity even
correction, the most general odd correction can be obtained by replacing
$\mathcal{G}$ with:

$\widetilde{\mathcal{G}}:=(\star\mathcal{C})_{abcd}\mathcal{C}^{abcd}\,.$
(117)

As argued at the end of section II.6, $\mathcal{P}$ anticommutes with
$\widetilde{\mathcal{G}}$:

$\mathcal{P}\widetilde{\mathcal{G}}+\widetilde{\mathcal{G}}\mathcal{P}\,=0,$
(118)

thus, the action of $\mathcal{P}$ in the EOM will be:

$\mathcal{P}\left(\mathcal{O}_{0}-L^{2}\lambda\widetilde{\mathcal{G}}\right)\Phi=\left(\mathcal{O}_{0}+L^{2}\lambda\widetilde{\mathcal{G}}\right)\mathcal{P}\Phi\,.$
(119)

$\mathcal{P}$ flips the sign of $\lambda$. Thus, again, we can relate QNMs
whose frequency has positive real part with the ones with negative real part.
If $(m,\omega(\lambda))$ defines a QNM frequency, then so does
$(-m,-\bar{\omega}(-\lambda))$. Bearing this in mind, and following the
approach in section III.1, we could obtain the parity odd EFT corrections to
scalar QNMs numerically.

###  III.3 Results

Using the pseudo-spectral methods outlined in section A.3, we solved the
background and perturbation EOM. In figure 4 we show that, as expected, these
methods converge exponentially with increasing grid size. We fixed the size of
the radial/angular grid and computed $\omega_{1}$ for varying angular/radial
grid sizes. By comparing the result with a more accurate value obtained for a
large grid, we see that the accuracy of $\omega_{1}$ increases exponentially.

The main results are presented in figure 1 and table 2. For the background
QNMs we obtain slightly unusual results. The imaginary part of the frequency
controls the decay rate of QNMs. The larger it is, the slower the decay rate.
Thus, the longest-lived QNM has the largest imaginary part of the frequency.
For gravitational and, as we will see below, electromagnetic QNMs of most BHs,
the longest-lived QNMs have $\ell=m=\absolutevalue{s}$. Usually, the ringdown
spectrum is mostly driven by these modes, with the relative intensity of modes
with smaller imaginary parts decaying exponentially. However, in figure 7, we
show that in the Schwarzschild limit the slowest decaying modes have
$\ell\rightarrow\infty$. This is consistent with the literature, (see table 1
of [37]). By inspecting figure 1, we show that this feature generalizes for
BHs with arbitrary spin. Note that this does not guarantee that the ringdown
does not include small $\ell$ modes. The non-linear dynamics of the BH
collision control the initial intensity of each QNM, and that could be heavily
skewed towards modes with small $\ell$. An in-depth discussion of this lies
beyond the scope of this paper.

In the Schwarzschild limit, the spacetime is spherically symmetric, and the
EOM can be solved non-perturbatively in $\widetilde{\lambda}$, using the
direct methods in section A.2. We obtained the QNM frequency $\omega$ for
several values of $\widetilde{\lambda}$. Taking $\widetilde{\lambda}$ to be
arbitrarily small, we disentangled $\omega$ into a background component
$\omega_{0}$ and the EFT correction $\omega_{1}$. The results can be seen in
the plot markers of figure 5. For $a\neq 0$, the EOM are no longer separable,
thus we proceed perturbatively, as outlined in section III.1. The results are
shown in the plot lines of figures 1 and 5, and in table 2. We found excellent
agreement between the direct results in the Schwarzschild limit and the ones
for slowly rotating Kerr BHs. This happens both for the background and EFT
correction of the QNM frequencies.

Except very close to extremality, the real and imaginary parts of $\omega_{1}$
have maximal absolute value for $\ell=m$ modes. For $\ell\neq 0$, the
correction is largest at $a\sim 0.9M$, quickly approaching $0$ afterwards. For
$a\gtrsim 0.9M$, the correction to modes with $\ell\neq m>0$, or $\ell=m=0$
increases quickly, before approaching finite values at $a\rightarrow M$. This
can more easily be seen in a parametric form in figure 5. $\ell=m\neq 0$
modes, also known as slowly damped QNMs, are particularly interesting in the
near extremal limit. Thus, in figure 6, we show the QNM frequencies in the
region $0.99<a/M\lesssim 0.999975$. As predicted in [38], we find that
$\omega_{0}$ approaches $\ell/2$. Furthermore, for $\ell=m=2$, the convergence
rate of $M\omega_{0}-1$ seems to be consistent with $\sqrt{1-a/M}$, as
predicted in [38]. However, for $\ell=m=1$, $M\omega_{0}$ seems to approach to
$1/2$ slightly faster than expected. The author is convinced that the
behaviour should agree better with the prediction in [38] if we studied BHs
even closer to extremality. However, extracting accurate results in this
regime involves using very large collocation grids, increasing dramatically
the size of the linear systems to solve. Thus, we did not cross $a\sim
0.999975M$. On the other hand, we found that, for slowly damped modes,
$M\omega_{1}$ approaches $0$ in the same limit. As outlined in the caption of
figure 6, the near-extremal region of parameter space has an intricate
structure. Both the real and imaginary parts of $\omega_{1}$ cross $0$ very
close to extremality, making it very difficult to accurately estimate the
convergence rate. Again, this might be overcome by studying BHs that are even
closer to extremality.

##  IV Corrections to electromagnetic modes

Drawing from the methods and intuition developed in section III, we will
obtain the corrections to electromagnetic QNMs, under the parity even / odd
terms introduced in the action (6). Contrary to a scalar field,
electromagnetism has two degrees of freedom, encoded in the two polarizations
of light. In a vacuum, both follow null geodesics of the background spacetime,
propagating at the speed of light. However, when passing through certain
media, the propagation speed can be polarization dependent, leading to
birefringent behaviour. In [50], considering the parity even part of 6, the
authors proved that the EFT correction, effectively turns vacuum into a
birefringent medium. For a type-D spacetime, in the Eikonal limit, they found
that light propagates along the null geodesics of two different effective
metrics, depending on its polarization. There is a known correspondence
between null geodesics and the large $\ell$ limit of QNM frequencies, see [13,
14, 15]. Hence, the birefringent behaviour of light should lead to a splitting
of the QNM frequencies into two families. In section IV.2, we will set
$\lambda_{(o)}=0$ in the action (6), and obtain the corrections to the two
families under parity even EFT corrections. Then, in section IV.3, we will
keep $\lambda_{(o)}\neq 0$ and generalize the result to mixed parity EFT
corrections. Surprisingly, we find that the correction to QNM frequencies due
to parity odd EFT terms coincides with the parity even case. We find that, for
each background QNM frequency there are two possible corresponding corrections
with opposing signs:

$\omega^{\pm}=\omega_{0}+\left(\frac{L}{M}\right)^{2}\absolutevalue{\lambda}\,\omega_{1}^{\pm}\,.$
(120)

where:

$\omega_{1}^{\pm}=\pm(-1)^{\ell+m+1}\,\omega_{1}$ (121)

and:

$\absolutevalue{\lambda}=\sqrt{\lambda_{(e)}^{2}+\lambda_{(o)}^{2}}\,.$ (122)

The factor of $(-1)^{\ell+m+1}$ was introduced to ensure that, in all cases
studied, the Schwarzschild limit of $\omega_{1}$ has positive real and
imaginary parts.666In section IV.2.1, we will show that, for parity even
corrections, in the Schwarzschild limit, the imaginary part of $\omega_{1}$
vanishes in the limit $\ell\rightarrow\infty$. As will be clear in section
IV.4.1, this feature does not generalize for finite $\ell$, thus we can, in
general, enforce the positivity of $\mathcal{I}m\,\omega_{1}$. Because
$\omega_{1}^{+}=-\omega_{1}^{-}$, this sign carries no physical meaning. This
is in close analogy to the results for the gravitational case, see [33, 34,
35]. Finally, in section IV.4, we will use the pseudo-spectral methods in the
appendix A to solve the EOM, obtaining the corrections to the QNM frequencies
in all cases.

###  IV.1 Equations of motion

In this section, we will derive the corrected EOM for electromagnetic waves in
a Kerr background. Using $p$ forms and hodge duals, we will cast them into a
form that can be easily tackled in later sections. Varying action (6) with
respect to the vector potential, and using equation (81), we get:

$\nabla^{b}\left(F_{ab}-\frac{1}{3}L^{2}\left(\lambda_{(e)}\mathcal{C}+\lambda_{(o)}\star\mathcal{C}\right)_{ab}{}^{cd}F_{cd}\right)=0\,.$
(123)

We want to solve this equation directly in terms of the components of
$F_{ab}$, thus, we must also enforce the electromagnetic Bianchi identity,
$\differential F=0$, guaranteeing $F$ is a closed form. Following the approach
in section II.7, the two conditions can be encoded by the real and imaginary
part of a single complex equation. We get:

$i\star\differential\left(F^{-}+i\frac{1}{3}L^{2}\left(\lambda_{(e)}\star\mathcal{C}-\lambda_{(o)}\mathcal{C}\right)\,:F\right)=0\,,$
(124)

where, $F^{-}$ was defined in equation (82), and we define “$:$” as a double
contraction. For some 4 tensor $\mathcal{A}_{ab}{}^{cd}$ and 2 tensor
$\mathcal{B}_{ab}$ we have:

$\left(\mathcal{A}:\mathcal{B}\right)_{ab}:=\mathcal{A}_{ab}{}^{cd}\mathcal{B}_{cd}\,.$
(125)

Assuming $\mathcal{C}$ and $F$ are real, the real part of (124) reduces to
equation (123), whereas the imaginary part encodes the Bianchi identity. By
solving the equation for a general $F$, we guarantee that both conditions are
satisfied. This equation takes a simpler form, if we aggregate $\lambda_{(e)}$
and $\lambda_{(o)}$ into a single complex coupling constant:

$\lambda=\lambda_{(e)}+i\lambda_{(o)}\,.$ (126)

Using the identities in section II.6, we get:

$i\star\differential\left(F^{-}-L^{2}\frac{\bar{\lambda}}{6}\,\mathcal{C}^{-}:F^{-}+L^{2}\frac{\lambda}{6}\,\mathcal{C}^{+}:F^{+}\right)=0\,.$
(127)

By definition, $F^{+}=\bar{F}^{-}$, thus the EOM depend on both, $F^{-}$ and
its complex conjugate. Consequently, the equation is not holomorphic in
$F^{\pm}$, and we must independently solve its real and imaginary parts.

To make progress, define $\mathcal{F}^{\pm}$, such that:

$\mathcal{F}^{\pm}=F^{\pm}-\left(\frac{L}{M}\right)^{2}\left(\lambda_{(e)}\pm
i\lambda_{(o)}\right)M^{2}\mathcal{C}^{\pm}:F^{\pm}\,.$ (128)

Note that $\mathcal{F}^{\pm}$ obeys the same properties as $F^{\pm}$. We have:

$\displaystyle\bar{\mathcal{F}}^{\pm}$ $\displaystyle=\mathcal{F}^{\mp}\,,$
(129) $\displaystyle\star\mathcal{F}^{\pm}$ $\displaystyle=\pm
i\mathcal{F}^{\pm}\,.$

Thus, we can define:

$\mathcal{F}=\mathcal{F}^{+}+\mathcal{F}^{-}\,,$ (130)

as a $2$ form that reduces to $F$ in the $L\rightarrow 0$ limit. Working
perturbatively in $L/M$, in terms of $\mathcal{F}^{\pm}$, the EOM are:

$\star\differential\left(\mathcal{F}^{-}+L^{2}\frac{\lambda}{6}\,\mathcal{C}^{+}:\mathcal{F}^{+}\right)=\mathcal{O}\left(\frac{L}{M}\right)^{4}\,.$
(131)

To project this into a null tetrad basis, we must first define NP scalars to
encode $\mathcal{F}$. Following the convention in equation (37), we get:

$\displaystyle\Phi_{0}$
$\displaystyle:=l^{a}m^{b}\mathcal{F}_{ab}=\left(1-\frac{1}{3}\bar{\lambda}\,L^{2}\,\psi_{2}\right)\phi_{0}\,,$
(132) $\displaystyle\Phi_{1}$
$\displaystyle:=\frac{1}{2}\left(l^{a}n^{b}-m^{a}\bar{m}^{b}\right)\mathcal{F}_{ab}=\left(1+\frac{2}{3}\bar{\lambda}\,L^{2}\,\psi_{2}\right)\phi_{1}\,,$
$\displaystyle\Phi_{2}$
$\displaystyle:=\bar{m}^{a}n^{b}\mathcal{F}_{ab}=\left(1-\frac{1}{3}\bar{\lambda}\,L^{2}\,\psi_{2}\right)\phi_{2}\,.$

In the limit $L\rightarrow 0$, $\Phi_{i}\rightarrow\phi_{i}$. The Weyl tensor
also obeys a Bianchi identity $\nabla_{[e}\mathcal{C}_{ab]cd}=0$. In a vacuum,
this implies that $\nabla^{a}\mathcal{C}_{abcd}=0$. Putting the two identities
together, we can prove:

$\left(\star\differential\,\mathcal{C}^{+}:\mathcal{F}^{+}\right)_{a}=\mathcal{C}^{+}_{ab}{}^{cd}\nabla^{b}\mathcal{F}_{cd}^{+}\,.$
(133)

Using this, we can prove that in the null tetrad basis, (131) reduces to:

$\mathcal{M}=\frac{1}{3}\bar{\psi}_{2}L^{2}\lambda\,J+\mathcal{O}\left(\frac{L}{M}\right)^{4}\,,$
(134)

where:

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

and:

$J_{(\alpha)}:=\begin{pmatrix}\left(\eth+2\bar{\tau}^{\prime}\right)\bar{\Phi}_{0}+2\left(\text{\thorn}+\bar{\rho}\right)\bar{\Phi}_{1}\\\
-\left(\eth^{\prime}+2\bar{\tau}\right)\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\left(\eth+\bar{\tau}^{\prime}\right)\bar{\Phi}_{1}\\\
\left(\text{\thorn}^{\prime}+2\bar{\rho}^{\prime}\right)\bar{\Phi}_{0}+2\left(\eth^{\prime}+\bar{\tau}\right)\bar{\Phi}_{1}\end{pmatrix}_{(\alpha)}\,.$
(136)

As expected from the identities $n^{a}=\left(l^{a}\right)^{\prime}$ and
$\bar{m}^{a}=\left(m^{a}\right)^{\prime}$, the $2^{nd}$ and $4^{th}$ line can
be obtained from the $1^{st}$ and $3^{rd}$ by using the prime operation.

While $M_{(\alpha)}$ depends exclusively on $\Phi_{i}$, $J_{(\alpha)}$, is a
function of $\bar{\Phi}_{i}$. Thus, if we try to factor out $\exp(-i\omega t)$
from each $\Phi_{i}$, we will get a term proportional to $\exp(-i\omega t)$ on
the LHS and a term proportional to $\exp(i\bar{\omega}t)$ on the RHS. The
standard approach to define QNMs cannot be applied directly. As remarked
before, this is a consequence of the non holomorphic character of equation
(127). Physically, this happens because the EFT correction breaks the
degeneracy between the two polarizations of light. To tackle this issue, we
must combine (134) with its complex conjugate, and choose a basis for
$\Phi_{i}$ and $\bar{\Phi}_{i}$, that leads to decoupled EOM. This is
analogous to finding a basis for the polarizations of light that propagate
along the two effective light cones of the corrected theory. If we consider
exclusively parity even corrections, we expect that this basis will be related
with parity eigenstates. In section IV.2, we will make this statement more
concrete. When considering mixed parity corrections, the answer is not as
obvious. We will explore this case in section IV.3.

###  IV.2 Parity even corrections

In this section, we set $\lambda_{(o)}=0$ and consider only parity even EFT
corrections. In section IV.2.1, Using the effective metrics derived in [22],
we will obtain an analytic expression for the QNMs in the large $\ell$ limit.
Then, in section IV.2.2, we specialize the discussion to the Schwarzschild
limit. Working non-perturbatively in $L/M$, we obtain decoupled ODEs for the
QNMs that can readily be solved with the direct methods in the appendix A. We
compare these equations with the ones derived in [51], and show that they are
equivalent. Finally, in section IV.2.3, we generalize the discussion to Kerr
BHs of any spin, obtaining decoupled PDEs for the QNMs, that can be solved
with the pseudo-spectral approach outlined in A.3.

####  IV.2.1 Eikonal limit

The Eikonal approximation, also known as the geometric optics limit, is valid
when the wavelength of light, dubbed $\varepsilon$, is much smaller than the
characteristic length scale of the background spacetime. For a Kerr
background, this is the statement that $\alpha\ll 1$, with:

$\alpha=\frac{\varepsilon}{M}\,.$ (137)

In [22], the authors explored this limit for the parity even part of (6),
considering several algebraically special BH spacetimes. In the beginning of
this section, we will outline their approach, and obtain the effective metrics
that encode the propagation of light on a Kerr background. Then, we will
exploit the connection with the large $\ell$ limit of QNMs to derive analytic
expressions for the frequencies in this limit.

We start with a plane wave Ansatz:

$A^{a}=\mathcal{A}^{a}e^{i\frac{\Theta}{\alpha}}\,.$ (138)

We assume that the characteristic length scale of all variables is $M$, with
the factor of $1/\alpha$ in the exponential ensuring that the signal has a
short wavelength. Plugging this into equation (123), and setting
$\lambda_{(o)}$ we get:

$C_{abcd}K^{b}K^{c}\mathcal{A}^{d}+\mathcal{O}\left(\frac{\varepsilon}{M}\right)=0\,.$
(139)

where:

$C_{ab}{}^{cd}=\delta_{a}^{[c}\,\delta_{b}^{d]}-\frac{1}{3}\,\lambda_{(e)}L^{2}\,\mathcal{C}_{ab}{}^{cd}\,,$
(140)

$K_{a}:=\nabla_{a}\Theta$ is the wavevector of (138), and parametrizes the
propagation direction.

As pointed out in [66], there is an issue with this approach. By dimensional
analysis, we have $\mathcal{C}\sim 1/M^{2}$, thud, the size of the correction
in (140) is:

$\left(\frac{L}{M}\right)^{2}\lambda_{(e)}=\left(\frac{L}{\varepsilon}\right)^{2}\left(\frac{\varepsilon}{M}\right)^{2}\lambda_{(e)}\,.$
(141)

Now, the EFT approximation is only valid, if $L$ is much smaller than any
other length scale in the theory. In particular, $L/\varepsilon\ll 1$, or
else, higher order derivative corrections become relevant. Thus:

$\left(\frac{L}{M}\right)^{2}\ll\left(\frac{\varepsilon}{M}\right)^{2}$ (142)

In the Eikonal approximation, we neglected all terms that are
$\mathcal{O}(\varepsilon/M)$, thus we must also neglect the contribution of
the EFT correction. The EFT correction is invisible in the Eikonal limit.

The same argument can be made on the QNM side. Bellow, we will obtain an
asymptotic expansion of the QNM frequencies in the large $\ell$ limit.
Roughly, we have:

$\omega=\ell\,\omega_{(1)}+\omega_{(0)}+\frac{1}{\ell}\omega_{(-1)}+\cdots$
(143)

The angular part of a QNM with azimuthal number $\ell$ has $\ell$ nodes, thus,
the length-scale of this mode is at least $M/\ell$. Hence, the validity of the
EFT approximation requires $L/M\ll 1/\ell$. We expect the correction to the
QNM frequency to be $\mathcal{O}(L/M)^{2}$ smaller than the background, thus
it should only be present at $\mathcal{O}(1/\ell^{2})$ above. We conclude that
we cannot extract information from the Eikonal limit while being consistent
with the EFT approximation. In spite of these facts, throughout the remainder
of this section, we compute the EFT corrections to large $\ell$ QNMs, ignoring
the requirement $L\ll\varepsilon$. We will do so exclusively to derive an
analytical check to test our numerics.

$C_{ab}{}^{cd}$, usually known as the susceptibility tensor, is a rank-4
tensor with the same algebraic symmetries as the Weyl tensor. Thus, equation
(139) is trivially solved if $\mathcal{A}^{a}$ is in the same direction as
$K^{a}$, however, this is a pure gauge solution
($A_{a}\sim\nabla_{a}e^{i\Theta}$), and thus, unphysical. To obtain physical
solutions, we must pick the direction of $K$ that the matrix
$U_{ad}=C_{abcd}K^{b}K^{c}$ has at most rank 2. This requirement leads to the
derivation of the Fresnel equation. This is a quartic equation on $K$, that
takes the form:

$\mathcal{G}_{abcd}K^{a}K^{b}K^{c}K^{d}=0\,.$ (144)

where $\mathcal{G}^{abcd}$ is a cubic contraction of $C_{abcd}$ defined in
equation 15 of [50].

For a type-D spacetime, choosing a frame aligned with the two principal null
directions, The Fresnel equation factorizes as the product of two equations,
quadratic in $K$. We get [50]:

$\left(\mathfrak{g}^{(+)}_{ab}K^{a}K^{b}\right)\left(\mathfrak{g}^{(-)}_{cd}K^{c}K^{d}\right)=0\,$
(145)

where:

$\mathfrak{g}^{(\pm)}_{ab}=g_{ab}+\sigma^{\pm}\left(m_{a}\bar{m}_{b}+m_{b}\bar{m}_{a}\right)\,.$
(146)

Here, $\mathfrak{g}^{(\pm)}_{ab}$ is the background Kerr metric, and:

$\sigma^{\pm}=\frac{\pm
2\lambda_{(e)}L^{2}\left|\psi_{2}\right|}{\left|1+\frac{1}{3}\lambda_{(e)}L^{2}\left(2\psi_{2}-\bar{\psi}_{2}\right)\right|\mp\lambda_{(e)}L^{2}\left|\psi_{2}\right|}\,.$
(147)

For $L\ll M$, $\mathfrak{g}^{(\pm)}_{ab}$ are non-degenerate symmetric
2-tensors, with Lorentzian signature. Thus, they can be understood as
effective metrics, with the corresponding solutions $K^{a}$ following their
null curves. In fact, because $K_{a}$ is the gradient of a scalar, we can
prove that these curves must be geodesic. Thus, in the Eikonal limit,
depending on the polarization, light will propagate along either of two
distinct effective light cones. The EFT correction effectively turns spacetime
into a birefringent medium.

To draw a parallel between the analytic results obtained here and the
numerical results obtained in section IV.4.2, we will work perturbatively in
$L/M$. At leading order, we have:

$\sigma^{\pm}=\pm
2\lambda_{(e)}L^{2}\absolutevalue{\psi_{2}}+\mathcal{O}\left(\frac{L}{M}\right)^{4}\,.$
(148)

Substituting this in equation (146), we get:

$\mathfrak{g}^{(\pm)}_{ab}=g_{ab}+2\widetilde{\lambda}_{(\pm)}M^{2}\absolutevalue{\psi_{2}}\,\left(m_{a}\bar{m}_{b}+m_{b}\bar{m}_{a}\right)\,,$
(149)

where:

$\widetilde{\lambda}_{(\pm)}:=\pm\left(\frac{L}{M}\right)^{2}\lambda_{(e)}\,.$
(150)

For a Kerr spacetime, the large $\ell$ limit of QNMs is fully specified by
circular photon orbits [13, 14, 15]. In fact, $\ell=\pm m$ QNMs are related to
the two sets of unstable equatorial null geodesics. The following relation
holds:

$\omega_{\ell=\pm
m}=\ell\,\Omega+\mathcal{O}\left(\ell^{0}\right)+i\left(\left(\frac{1}{2}+n\right)\gamma^{L}+\mathcal{O}\left(\ell^{-2}\right)\right)\,.$
(151)

Here, $\Omega=\frac{d\phi}{dt}$ represents the angular frequency of the orbit,
and $\gamma^{L}$ is the corresponding Lyapunov exponent, which indicates the
instability time scale with respect to the coordinate $t$. The order of
magnitude of the corrections was derived in [14], by comparing the Eikonal
limit with a WKB expansion of the Teukolsky equation. We will make the logical
leap that the results generalize to photons propagating in the effective
metrics (146). Intuitively, this should be the case because of the
similarities between the geometric optics approximation and a WKB expansion,
however, without a more formal treatment, we cannot make a statement regarding
the quality of the approximation. $\ell=m$ modes relate to null-geodesics that
co-rotate with the BH, while $\ell=-m$ modes are specified by counter-rotating
photons.

(a)

(b)

Figure 8: Here, we show the prediction for the large $\ell$ limit of
electromagnetic QNM frequencies. On the left, we have the background
component, and on the right, we present the corresponding correction, as
defined in equation (120). For each plot, the left side showcases the
$\ell=-m$ modes, with $a/M\in[0,1]$, running from the middle to the left.
Conversely, the right side illustrates the $\ell=m$ modes, with $a/M\in[0,1]$,
from the middle to the right. We arrange the plots in this manner to
underscore the equivalence between the $m\rightarrow-m$ operation and
$a\rightarrow-a$. In the near extremal limit, for $\ell=m$ modes, we find that
$M\omega_{0}^{\ell=m}\rightarrow\ell/2$, and $M\omega_{1}^{\ell=m}\rightarrow
0$. Finally, note that in the Schwarzschild limit,
$\mathcal{I}m(\omega_{1})\rightarrow 0$, as can be deduced from equation (165)

Perturbatively, the two effective metrics are related by reversing the sign of
$\lambda$. Thus, in the Eikonal limit, the corresponding QNMs must be related
by reversing the sign of $\lambda$. Thus, following the convention in equation
(120), we can conclude that:

$\omega^{\absolutevalue{\ell}=\absolutevalue{m}\rightarrow\infty}_{\pm}=\omega_{0}\mp\left(\frac{L}{M}\right)^{2}\lambda_{(e)}\,\omega_{1}+\mathcal{O}\left(\frac{L}{M}\right)^{4}\,.$
(152)

Null geodesics extremize the Lagrangian,

$\mathcal{L}=\frac{1}{2}\mathfrak{g}^{(\pm)}_{ab}\derivative{x^{a}}{\chi}\derivative{x^{b}}{\chi}\,.$
(153)

where $\chi$ is some affine parameter. By definition $\mathcal{L}=0$,
furthermore due to the symmetries of $\mathfrak{g}$, we have that:

$\displaystyle\mathcal{E}=\partialderivative{\mathcal{L}}{\left(\differential
t/\differential\chi\right)},$ (154) $\displaystyle
L=\partialderivative{\mathcal{L}}{\left(\differential\phi/\differential\chi\right)}\,,$

are conserved quantities. Inverting this relation, we can express
$(\differential t/\differential\chi,\differential\phi/\differential\chi)$ in
terms of $(\mathcal{E},L)$. Because we are looking for equatorial geodesics,
$x=\differential x/\differential\chi=0$. Substituting all of this in
$\mathcal{L}=0$, and solving with respect to $\differential
r/\differential\chi$, we get:

$\left(\frac{1}{\mathcal{E}}\derivative{r}{\chi}\right)^{2}=V(r,b):=V_{0}(r,b)+\widetilde{\lambda}_{(\pm)}V_{1}(r,b)\,,$
(155)

where:

$\displaystyle V_{0}(r,b)$
$\displaystyle=\frac{2(a+b)^{2}\,M+(a^{2}-b^{2})\,r+r^{3}}{4r^{3}}\,,$ (156)
$\displaystyle V_{1}(r,b)$
$\displaystyle=\frac{(a+b)^{2}\Delta_{r}}{2r^{7}}M^{2}\,,$

and:

$b=\frac{L}{\mathcal{E}}\,$ (157)

is the impact parameter of the null geodesic. Taking the derivative of (155)
with respect to $\chi$, we see that circular geodesics happen if and only if:

$\begin{cases}V(r,b)=0\,,\\\ \partialderivative{r}V(r,b)=0\,.\end{cases}$
(158)

To solve this, we will make the Ansatz:

$\displaystyle r$
$\displaystyle=r_{0}+\widetilde{\lambda}_{(\pm)}r_{1}+\mathcal{O}\left(\widetilde{\lambda}_{(\pm)}\right)^{2}\,,$
(159) $\displaystyle b$
$\displaystyle=b_{0}+\widetilde{\lambda}_{(\pm)}b_{1}+\mathcal{O}\left(\widetilde{\lambda}_{(\pm)}\right)^{2}\,.$

At $\mathcal{O}\left(\widetilde{\lambda}_{(\pm)}{}^{0}\right)$, the equations
reduce to:

$\begin{cases}r_{0}\,(r_{0}-3M)^{2}=4a^{2}M\,,\\\
b_{0}=a\frac{r_{0}+3M}{r_{0}-3M}\,.\end{cases}$ (160)

The cubic equation has two roots for $r>r_{+}$, given by [14]:

$r_{0}^{\left(\pm\right)}=2M\left(1+\cos\left(\frac{2}{3}\cos^{-1}\left(\mp\frac{a}{M}\right)\right)\right)$
(161)

Notice that $r_{0}^{+}\leq 3M\leq r_{0}^{-}$, with equality only for $a=0$.
The $+$ case describes photons that co-rotate with the BH, being dual to QNMs
with positive real part, whereas the $-$ case describes counter-rotating
photons, dual to QNMs with negative real part. This $\pm$ is independent from
the one in $\widetilde{\lambda}_{(\pm)}$, thus, for clarity, we will omit it
below.

The first order is a linear equation for $(r_{1},b_{1})$, thus we can solve it
directly. Using equation (160) to simplify the resulting expression, we get:

$\displaystyle b_{1}$
$\displaystyle=M\,\frac{(r_{0}-3M)(r_{0}-M)}{4\,a\,r_{0}}\,,$ (162)
$\displaystyle r_{1}$
$\displaystyle=M^{2}\,\frac{5a^{2}\,+\,r_{0}\,(4\,r_{0}-9\,M)}{3\,r_{0}^{3}}\,.$

We can now compute the relevant parameters for equation (151). The Lyapunov
exponent is defined as:

$\gamma^{L}:=\sqrt{\frac{\mathcal{E}^{2}\,\partial^{2}V(r,b)/\partial
r^{2}}{(\differential t/\differential\chi)^{2}}}\,.$ (163)

Expanding in powers of $\widetilde{\lambda}_{(\pm)}$, we get:

$\gamma^{L}=\gamma^{L}_{0}+\widetilde{\lambda}_{(\pm)}\gamma^{L}_{1}+\mathcal{O}\left(\widetilde{\lambda}_{(\pm)}\right)^{2}$
(164)

with:

$\displaystyle\gamma^{L}_{0}$
$\displaystyle=\frac{\sqrt{3}}{r_{0}}\frac{r_{0}-M}{r_{0}+3\,M}\,,$ (165)
$\displaystyle\gamma^{L}_{1}$
$\displaystyle=a\,\gamma^{L}_{0}M\frac{(r_{0}-M)(5\,r_{0}+9\,M)}{12\,b_{0}r_{0}^{3}}\,.$

Similarly, the orbital frequency is given by:

$\Omega:=\frac{\differential\phi/\differential\chi}{\differential
t/\differential\chi}=-\frac{1}{b_{0}}\,+\,\widetilde{\lambda}_{(\pm)}\frac{\gamma^{L}_{0}}{\sqrt{3}}{\,b_{0}\,r_{0}}+\mathcal{O}\left(\widetilde{\lambda}_{(\pm)}\right)^{2}$
(166)

Plugging these results into equation (151), and using the convention in
equation (152), we get an analytic estimate for the correction to QNM
frequencies in the large $\ell$ limit. The results are summarized in figure 8.

####  IV.2.2 Schwarzschild limit

In the Schwarzschild limit, considering parity even EFT corrections, decoupled
EOM for the QNMs were first obtained in [51]. There, the authors worked
directly with the vector potential $A_{a}$ and obtaining the EOM using vector
spherical harmonics. In this section, we will arrive at the same equations
independently, by working with the gauge invariant components of $F_{ab}$,
using the GHP formalism. This way, it will be easier to relate our results
with the ones obtained later for general Kerr BHs in terms of the same
variables. For non-spinning BHs, the spacetime is spherically symmetric, thus,
the EOM are much simpler. This allows us to proceed non-perturbatively in
$L/M$, ignoring the EFT validity conditions. We will do so to relate our
results with the ones in [51], but we must bear in mind that only solutions
with $L\ll M$ are relevant from an EFT perspective. Spherical symmetry
simplifies the GHP formalism. In addition to the GHP identities that are
satisfied for any type-D spacetime (see equation (39)), we have:

$\displaystyle\rho=\bar{\rho},\quad\quad$
$\displaystyle\rho^{\prime}=\bar{\rho}^{\prime}$ (167)
$\displaystyle\psi_{2}=\bar{\psi}_{2},$ $\displaystyle\tau=\tau^{\prime}=0$

Furthermore, all scalars depend exclusively on the radial coordinate $r$.
Projecting equation (127) into a tetrad basis, and setting $\lambda_{o}=0$, we
get:

$\begin{pmatrix}\eth^{\prime}\phi_{0}-\left(\text{\thorn}-2\rho\right)\phi_{1}\\\
\eth\phi_{2}-\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)\phi_{1}\\\
\left(\text{\thorn}^{\prime}-\rho^{\prime}\right)\phi_{0}-\eth\phi_{1}\\\
\left(\text{\thorn}-\rho\right)\phi_{2}-\eth^{\prime}\phi_{1}\end{pmatrix}=\frac{1}{3}\lambda_{(e)}L^{2}\psi_{2}\begin{pmatrix}\eth^{\prime}\phi_{0}+\eth\bar{\phi}_{0}+2\left(\text{\thorn}+\rho\right)\left(\phi_{1}+\bar{\phi}_{1}\right)\\\
\eth\phi_{2}+\eth^{\prime}\bar{\phi}_{2}+2\left(\text{\thorn}^{\prime}+\rho^{\prime}\right)\left(\phi_{1}+\bar{\phi}_{1}\right)\\\
\left(\text{\thorn}^{\prime}+2\rho^{\prime}\right)\phi_{0}-\left(\text{\thorn}+2\rho\right)\bar{\phi}_{2}+2\,\eth\left(\phi_{1}-\bar{\phi}_{1}\right)\\\
\left(\text{\thorn}+2\rho\right)\phi_{2}-\left(\text{\thorn}^{\prime}+2\rho^{\prime}\right)\bar{\phi}_{0}+2\,\eth^{\prime}\left(\phi_{1}-\bar{\phi}_{1}\right)\end{pmatrix}$
(168)

To simplify these equations, our strategy is to express them as a system of
real equations of real variables. $\phi_{1}$ has spin weight $0$ (see section
II.2), thus its GHP weights are invariant under complex conjugation. That is
why, in the RHS of (168), we see terms that explicitly depend on the real and
imaginary part of $\phi_{1}$, without breaking GHP covariance. This would not
be possible for $\phi_{0}$ and $\phi_{2}$, as they have spin weight $1$ and
$-1$ respectively. To represent them with real, GHP-covariant variables, we
must first raise / lower their spin weight by acting with $\eth$ and
$\eth^{\prime}$ respectively. Define the polarized Maxwell GHP scalars as:

$\displaystyle P_{0}^{+}=\mathcal{R}e\left(\eth^{\prime}\phi_{0}\right),\,$
$\displaystyle P_{1}^{+}=\mathcal{R}e\left(\phi_{1}\right),\,$ $\displaystyle
P_{2}^{+}=\mathcal{R}e\left(\eth\phi_{2}\right)\,,$ (169) $\displaystyle
P_{0}^{-}=\mathcal{I}m\left(\eth^{\prime}\phi_{0}\right),\,$ $\displaystyle
P_{1}^{-}=\mathcal{I}m\left(\phi_{1}\right),\,$ $\displaystyle
P_{2}^{-}=\mathcal{I}m\left(\eth\phi_{2}\right)\,.$

By definition, for any GHP gauge choice, these variables are real. On a
similar note, the first two equations in (168) have spin-weight $0$, whereas
the others have spin-weight $1$ and $-1$ respectively. Thus, after applying
$\eth^{\prime}$ and $\eth$ to the $3^{rd}$ and $4^{th}$ equations, the real
and imaginary parts of the resulting system are GHP-covariant. Remarkably,
after extensive use of GHP identities [53], we can prove that the real part of
equation (168) depends exclusively on $P_{i}^{+}$. Similarly, the imaginary
part depends exclusively on $P_{i}^{-}$. Thus, the EOM for the two
polarizations decouple:

$\begin{pmatrix}\left(\text{\thorn}-2\rho\right)P_{1}^{+}-P_{0}^{+}\\\
\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)P_{1}^{+}-P_{2}^{+}\\\
\eth\,\eth^{\prime}\,P_{1}^{+}-\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)P_{0}^{+}\\\
\eth\,\eth^{\prime}\,P_{1}^{+}-\left(\text{\thorn}-2\rho\right)P_{2}^{+}\end{pmatrix}=\frac{1}{3}\lambda_{(e)}L^{2}\psi_{2}\begin{pmatrix}4\left(\text{\thorn}+\rho\right)P_{1}^{+}+2P_{0}^{+}\\\
4\left(\text{\thorn}^{\prime}+\rho^{\prime}\right)P_{1}^{+}+2P_{2}^{+}\\\
\left(\text{\thorn}^{\prime}+\rho^{\prime}\right)P_{0}^{+}-\left(\text{\thorn}+\rho\right)P_{2}^{+}\\\
\left(\text{\thorn}+\rho\right)P_{2}^{+}-\left(\text{\thorn}^{\prime}+\rho^{\prime}\right)P_{0}^{+}\end{pmatrix}\,,$
(170)

and:

$\begin{pmatrix}\left(\text{\thorn}-2\rho\right)P_{1}^{-}-P_{0}^{-}\\\
\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)P_{1}^{-}-P_{2}^{-}\\\
\eth\,\eth^{\prime}\,P_{1}^{-}-\left(\text{\thorn}^{\prime}-2\rho^{\prime}\right)P_{0}^{-}\\\
\eth\,\eth^{\prime}\,P_{1}^{-}-\left(\text{\thorn}-2\rho\right)P_{2}^{-}\end{pmatrix}=\frac{1}{3}\lambda_{(e)}L^{2}\psi_{2}\begin{pmatrix}0\\\
0\\\
4\eth\,\eth^{\prime}\,P_{1}^{-}+\left(\text{\thorn}^{\prime}+\rho^{\prime}\right)P_{0}^{-}+\left(\text{\thorn}+\rho\right)P_{2}^{-}\\\
4\eth\,\eth^{\prime}\,P_{1}^{-}+\left(\text{\thorn}^{\prime}+\rho^{\prime}\right)P_{0}^{-}+\left(\text{\thorn}+\rho\right)P_{2}^{-}\\\
\end{pmatrix}\,.$ (171)

In the Schwarzschild limit, þ and $\text{\thorn}^{\prime}$ reduce to $(t,r)$
derivative operators, while $\eth$ and $\eth^{\prime}$ act exclusively on the
angle variables. Thus, the only angular dependence in equations (170) and
(171), is encoded in the $(\eth\,\eth^{\prime})$ operator. There is a known
relationship between $(\eth\,\eth^{\prime})$ and the spherical harmonics
$Y_{\ell m}$ [55]. We have:

$\eth\,\eth^{\prime}Y_{\ell m}=-\frac{\ell(\ell+1)}{2r^{2}}Y_{\ell m}\,.$
(172)

Given $\text{\thorn}Y_{\ell m}=\text{\thorn}^{\prime}Y_{\ell m}=0$, a
decomposition of the scalars in spherical harmonics removes all angular
dependence. Similarly, we can factor out time dependence by multiplying out
$\exp(-i\omega t)$.777Both $Y_{\ell m}$ and $\exp(-i\omega t)$ are complex
valued, however, we argued that $P_{i}^{\pm}$ must be real valued functions.
Because the EOM are now real, we can work with a complex-valued Ansatz and in
the end simply extract its real part. $P_{0}^{\pm}$ and $P_{2}^{\pm}$ have
boost weight $1$ and $-1$ respectively. Thus, to work with $\\{0,0\\}$
scalars, we will explicitly factor out $\rho$ and $\rho^{\prime}$
respectively. The final Ansatz will be:

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

Equations (170) and (171) reduce to two systems of four first-order ODEs for
the three radial variables, an overdetermined description. For both
polarizations, 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,
reducing the problem to a system of two ODEs with two variables. Finally, we
define $u^{\pm}(r)$ and $v^{\pm}(r)$, such that:

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

with

$\displaystyle f$ $\displaystyle:=1-\frac{2M}{r}\,,$ (175)
$\displaystyle\Xi(r)$
$\displaystyle:=1+\,\frac{2}{3}\,\widetilde{\lambda}_{(e)}\,\left(\frac{M}{r}\right)^{3}\,,$

where we defined:

$\widetilde{\lambda}_{(e)}:=\left(\frac{L}{M}\right)^{2}\lambda_{(e)}$ (176)

The EOM reduce to:

$\begin{cases}\derivative{r}v^{\pm}+\Xi^{\pm 1}\frac{i\omega}{f}\,u^{\pm}=0\\\
\derivative{r}u^{\pm}+\left(\Xi^{\mp
1}\frac{i\omega}{f}+\frac{1}{i\omega}V_{\pm}\right)v^{\pm}=0\end{cases}\,,$
(177)

where:

$\displaystyle V_{\pm}:=\frac{\ell(\ell+1)}{r^{2}}\zeta^{\mp 1}\,,$ (178)
$\displaystyle\zeta:=1-\frac{4}{3}\widetilde{\lambda}_{(e)}\left(\frac{M}{r}\right)^{3}\,.$

These equations can be combined into a single decoupled second-order ODE for
$v^{\pm}$. We get:

$\derivative{r}\Xi^{\mp 1}\,f\,\derivative{r}v^{\pm}\,+\,\left(\Xi^{\mp
1}\,\frac{\omega^{2}}{f}-V_{\pm}\right)v^{\pm}=0\,.$ (179)

After suitable variable changes, we can prove that this equation agrees with
the effective potentials derived in [51] (equations 16 and 17).888To match the
conventions in that paper, we must take $\lambda\rightarrow 12\,\alpha/L^{2}$.
In the $\widetilde{\lambda}_{(e)}\rightarrow 0$ limit, $\Xi\rightarrow 1$,
$\zeta\rightarrow 1$, and the EOM reduce to the wave equation for
electromagnetic fields in a Schwarzschild BH background, see [67].

As outlined in section II.4, to find the QNM spectrum of (179), we must set
BCs that are ingoing at $r\rightarrow r_{+}$, and outgoing at
$r\rightarrow\infty$. $\Xi$ and $\zeta$ are smooth at $r=r_{+}$, and asymptote
to $1+\mathcal{O}(1/r^{3})$ at $r\rightarrow\infty$. Thus, at leading order, a
Frobenius expansion near the boundaries is independent of
$\widetilde{\lambda}_{(e)}$. Hence, to set the appropriate BCs, we can proceed
exactly as in section II.4. Define $\widetilde{v}^{\pm}(r)$:

$\widetilde{v}^{\pm}(r)=\frac{1}{\alpha^{a=0}_{r}(r)}v^{\pm}(r)\,$ (180)

with:

$\alpha^{a=0}_{r}(r)=f^{-2i\omega M}\,\left(\frac{r}{M}\right)^{2i\omega
M}\,e^{ir\omega}$ (181)

We change variables to the compactified radial coordinates $y=1-r/r_{+}$, and
plug this in equation (179). QNMs are solutions to the resulting equation that
are smooth for $y\in[0,1]$. Treating $\widetilde{\lambda}_{(e)}$ as a free
parameter. this equation can be solved numerically, using pseudospectral
methods, as outlined in the appendix A. The EOM are not perturbative in
$\widetilde{\lambda}_{(e)}$, and can, in theory, be solved for
$\widetilde{\lambda}_{(e)}\in(-12,6)$.999This is the range such that $\zeta$
and $\Xi$ have no roots for $r>2M$. However, the EFT approximation entails
$\widetilde{\lambda}_{(e)}\ll 1$, thus, we should disregard the remaining of
the parameter space. Expanding equation (179) in powers of
$\widetilde{\lambda}_{(e)}$, we see that at leading order, the EOM for the
$\pm$ polarizations can be mapped into each other by changing
$\widetilde{\lambda}_{(e)}\rightarrow-\widetilde{\lambda}_{(e)}$. Thus, at
leading order, the two corrections to the QNM frequencies are equal and
opposite. This is in agreement with what we postulated in equation
(120).101010The convention we picked in this section actually differs from the
one in (120) by a factor of $(-1)^{(\ell+m)}$.

####  IV.2.3 Generalization to Kerr

In this section, we will solve equation (134), considering only parity even
EFT corrections (i.e. $\lambda_{(o)}=0$), for a Kerr BH background of any
spin. To do so, we will make explicit use of the EFT validity condition, and
work perturbatively in $L/M$, or equivalently, in terms of
$\widetilde{\lambda}_{(e)}$, defined in (176).

Take a generic differential equation:

$Q_{(0)}+\widetilde{\lambda}_{(e)}\left(\mathcal{G}Q_{(0)}+Q_{(1)}\right)+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})=0,$
(182)

where $Q_{(0)}$, $\mathcal{G}$ and $Q_{(1)}$ are arbitrary. At
$\mathcal{O}(\widetilde{\lambda}_{(e)}^{0})$, we have
$Q_{(0)}=\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})$. Substituting this on the
$\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})$ part of the equation, we get:

$Q_{(0)}+\widetilde{\lambda}_{(e)}Q_{(1)}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})=0,$
(183)

Thus, any identity that is satisfied at
$\mathcal{O}(\widetilde{\lambda}_{(e)}^{0})$ can be used to simplify the
$\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})$ terms. Following this idea, we
equate $\mathcal{M}_{(\alpha)}=0$ and solve with respect to the GHP
derivatives of $\phi_{1}$. The result can be used to eliminate all derivatives
of $\bar{\phi}_{1}$ from $J_{(\alpha)}$. Then, we apply the Teukolsky
$\zeta^{a}$ operator, defined in equation (91), to both sides of equation
(134). Using the GHP identities in [53] we get:

$\mathcal{O}_{1}\Phi_{0}=\lambda_{(e)}L^{2}\bar{\psi}_{2}\mathcal{J}_{i}\bar{\Phi}_{i}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\,,$
(184)

where the summation in $i$ is implicit, $\mathcal{J}_{i}$ are linear
operators:

$\displaystyle\mathcal{J}_{0}$
$\displaystyle=\left(2\left(\tau-2\bar{\tau}^{\prime}\right)-\eth\right)\eth+2\bar{\tau}^{\prime
2}\,,$ (185) $\displaystyle\mathcal{J}_{1}$
$\displaystyle=6\left(\rho\bar{\tau}^{\prime}+\bar{\rho}\left(\tau-4\bar{\tau}^{\prime}\right)\right)\,,$
$\displaystyle\mathcal{J}_{2}$
$\displaystyle=\left(2\left(\rho-2\bar{\rho}\right)-\text{\thorn}\right)\text{\thorn}+2\bar{\rho}^{2}\,,$

and $\mathcal{O}_{s}$ is the master Teukolsky operator defined in equation
(44).

From section II.7, we know that there is a Hertz potential $\psi_{H}$ such
that:

$\displaystyle\mathcal{O}_{-1}\psi_{H}=0\,,$ (186)
$\displaystyle\phi_{i}=\mathcal{T}_{i}\bar{\psi}_{H}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})\,,$

where $\mathcal{T}_{i}$ was defined in equation (104). Substituting this in
equation (184), we get:

$\begin{cases}\mathcal{O}_{1}\Phi_{0}=\widetilde{\lambda}_{(e)}\mathcal{D}\psi_{H}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\\\\[5.0pt]
\Phi_{0}=\mathcal{T}_{0}\bar{\psi}_{H}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})\end{cases}\,,$
(187)

with:

$\mathcal{D}:=M^{2}\bar{\psi}_{2}\,\mathcal{J}_{i}\,\bar{\mathcal{T}_{i}}\,.$
(188)

Equation (187) mixes modes proportional to $S$ with modes proportional to
$\bar{S}$. If we set $\Phi_{0}$ and $\psi_{H}$ proportional to
$S=\exp(-i\omega t+im\phi)$, $S$ factors out nicely from both sides of the
first equation in (187). However, for the second equation, we get $S$ on the
LHS and $\bar{S}$ on the RHS. The reverse problem occurs if we choose
$\psi_{H}$ to be proportional to $\bar{S}$ instead. In section II.5, we
studied the properties of the conjugate-parity-transform operator,
$\mathcal{P}$, defined in equation (66). We proved, that due to the parity
invariance of the Kerr metric, the operator commutes with the master Teukolsky
operator $\mathcal{O}_{s}$. Consequently, we found that if
$\psi^{(s)}=S(t,\phi)\Psi_{s}(r,x)$ is a QNM of $\mathcal{O}_{s}$ then so is
$\mathcal{P}\psi^{(s)}=(-1)^{m}\bar{S}(t,\phi)\Psi_{s}(r,-x)$. Taking linear
combinations of these states, we can express the kernel of $\mathcal{O}_{s}$
in a basis of $\mathcal{P}$ eigenstates. We restricted to parity invariant EFT
corrections, thus, we should be able to express solutions to equation (187) as
eigenstates of $\mathcal{P}$. Take:

$\displaystyle\Phi_{0}^{\pm}$
$\displaystyle:=\Phi_{0}\pm\mathcal{P}\Phi_{0}\,,$ (189)
$\displaystyle\psi_{H}^{\pm}$
$\displaystyle:=\psi_{H}\pm\mathcal{P}\phi_{H}\,.$

By definition, $\mathcal{P}\Phi_{0}^{\pm}=\pm\Phi_{0}^{\pm}$ and
$\mathcal{P}\psi_{H}^{\pm}=\pm\psi_{H}^{\pm}$, thus we associate these with
the polar and axial polarizations of light respectively.111111The polar
polarization is parity even, whereas its axial counterpart must be parity odd,
hence this classification. Because, $\mathcal{O}_{s}$, $\mathcal{D}$ and
$\mathcal{T}_{0}$ have even spin-weight, they all commute with $\mathcal{P}$:

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

Thus, applying $\mathcal{P}$ to both sides of equation (187), we get:

$\begin{cases}\mathcal{O}_{1}\,\mathcal{P}\Phi_{0}=\widetilde{\lambda}_{(e)}\mathcal{D}\,\mathcal{P}\psi_{H}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\\\\[5.0pt]
\mathcal{P}\Phi_{0}=\mathcal{T}_{0}\,\mathbb{P}\psi_{H}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})\end{cases}\,,$
(191)

where $\mathbb{P}$ is the parity transform operator, (see section II.5).
Writing equations (187) and (191) in the $\mathcal{P}$ eigenstates basis, we
get:

$\begin{cases}\mathcal{O}_{1}\Phi_{0}^{\pm}=\widetilde{\lambda}_{(e)}\mathcal{D}\psi_{H}^{\pm}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\\\\[5.0pt]
\Phi_{0}^{\pm}=\pm\mathcal{T}_{0}\,\mathbb{P}\psi_{H}^{\pm}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})\end{cases}\,.$
(192)

Because $[\mathcal{P},\mathcal{O}_{-1}]=0$, we have
$\mathcal{O}_{-1}\psi_{H}^{\pm}=0$. However, in a Kerr background, QNMs of the
Teukolsky equation are not, in general, eigenstates of $\mathcal{P}$,
consequently $\psi_{H}^{\pm}$ cannot be a QNM solution. The same can be said
about the $\mathcal{O}(\widetilde{\lambda}_{(e)}^{0})$ component of
$\Phi_{0}^{\pm}$. Nevertheless, we can choose $\psi_{H}^{\pm}$ and
$\Phi_{0}^{\pm}$ to be a linear combination of modes proportional to $S$ and
$\bar{S}$. Extending this to first order in $\widetilde{\lambda}_{(e)}$, we
take:

$\displaystyle\Phi_{0}^{\pm}$ $\displaystyle=\Phi\pm\mathcal{P}\Phi\,,$ (193)
$\displaystyle\psi_{H}^{\pm}$ $\displaystyle=\Psi\pm\mathcal{P}\Psi\,,$

with $\Psi$ and $\Phi$ proportional to $S$:121212It may seem that
$\Phi=\Phi_{0}$ and $\Psi=\psi_{H}$, however that is not the case. As argued
below equation (188), $(\Phi_{0},\psi^{H})$ cannot be proportional to $S$ or
$\bar{S}$, thus, necessarily $\Phi$ and $\Psi$ are different variables.

$\displaystyle\Psi(t,r,x,\phi)$ $\displaystyle=S(t,\phi)\,\Psi_{r,x}(r,x)\,,$
(194) $\displaystyle\Phi(t,r,x,\phi)$
$\displaystyle=S(t,\phi)\,\Phi_{r,x}(r,x)\,.$

$\mathcal{P}\Psi$ and $\mathcal{P}\Phi$ will be proportional to $\bar{S}$,
thus substituting this in equation (192) we get:

$\hskip-11.38092pt\begin{cases}S(t,\phi)\,\mathcal{F}^{\pm}(r,x)=\mp(-1)^{m}\bar{S}(t,\phi)\,\mathcal{P}\mathcal{F}^{\pm}(r,x)\\\\[5.0pt]
S(t,\phi)\,\mathcal{G}^{\pm}(r,x)=\mp(-1)^{m}\bar{S}(t,\phi)\,\mathcal{P}\mathcal{G}^{\pm}(r,x)\end{cases}\,,$
(195)

where:

$\displaystyle\mathcal{F}^{\pm}(r,x)$
$\displaystyle:=\mathcal{O}_{1,\omega}\Phi_{r,x}(r,x)-\widetilde{\lambda}_{(e)}\mathcal{D}_{\omega}\,\Psi_{r,x}(r,x)+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\,,$
(196) $\displaystyle\mathcal{G}^{\pm}(r,x)$
$\displaystyle:=\Phi_{r,x}\mp(-1)^{m}\,\mathcal{T}_{0,\omega}\,\mathbb{P}\Psi_{r,x}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})\,,$

with:

$\displaystyle\mathcal{D}_{\omega}$ $\displaystyle=\frac{1}{S}\mathcal{D}S\,,$
(197) $\displaystyle\mathcal{T}_{0,\omega}$
$\displaystyle=\frac{1}{S}\mathcal{T}_{0}S\,,$
$\displaystyle\mathcal{O}_{1,\omega}$
$\displaystyle=\frac{1}{S}\mathcal{O}_{1}S\,.$

$S$ and $\bar{S}$ are linearly independent, thus, their pre-factors must
vanish. We get:

$\displaystyle\mathcal{F}^{\pm}=0\,,$ (198a)
$\displaystyle\mathcal{G}^{\pm}=0\,.$ (198b)

Equation (198b) relates the background of $\Phi_{r,x}$ with the background of
$\Psi_{r,x}$. We choose the latter to be a QNM of the $s=-1$ Teukolsky
equation, and consequently a Hertz potential for the background QNM solutions.
This implies that, in the limit $\widetilde{\lambda}_{(e)}\rightarrow 0$,
$\Phi$ must be a QNM of the $s=1$ Teukolsky equation. We substitute this in
equation (198a), and work perturbatively in $\widetilde{\lambda}_{(e)}$ to
find the QNM spectrum of the full theory. Concretely, we will solve (198a) in
two steps. First, treating $\omega$ as a parameter independent of
$\widetilde{\lambda}_{(e)}$, we set the appropriate QNM BCs. As outlined in
section II.4, this is done by factoring out the appropriate corresponding
singular behaviour. Then, we expand $\omega$ in powers of
$\widetilde{\lambda}_{(e)}$, and solve equation (198b) perturbatively. If we
don’t follow this sequence, we cannot ensure that the BCs for QNMs are
properly enforced, and we might get logarithmic singularities plaguing the
numerics.

The singular behaviour of $\Psi_{r,x}$ is well understood. Working with the
compactified radial coordinate, we have:

$\Psi_{r,x}(r(y),x)=\alpha_{-1}(y,x)\widetilde{\Psi}_{y,x}(y,x)\,,$ (199)

where $\alpha_{-1}$, defined in (62), is singular at the domain border, and
$\widetilde{\Psi}_{y,x}$ is smooth for $(y,x)\in[0,1]\crossproduct[-1,1]$.
Working perturbatively in $\widetilde{\lambda}_{(e)}$, $\Phi_{r,x}$ takes the
form:

$\Phi_{r,x}(r(y),x)=\left(\alpha_{1}(y,x)+\widetilde{\lambda}_{(e)}\,\delta\alpha(y,x)\right)\left(\widetilde{\Phi}_{y,x}^{(0)}(y,x)+\widetilde{\lambda}_{(e)}\Phi_{y,x}^{(1)}(y,x)\right)+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\,.$
(200)

By definition, $\widetilde{\Phi}_{y,x}^{(0)}$ is smooth for
$(y,x)\in[0,1]\crossproduct[-1,1]$. For $\Phi$ to be a QNM, we require
$\widetilde{\Phi}_{y,x}^{(1)}$ to be smooth in the same interval.
$\delta\alpha$ is included to ensure that this is the case. It acts as a
counter-term to cancel non-smooth behaviour. With a Frobenius analysis of
equation (203), we can actually prove that $\delta\alpha=0$.131313This is
expected if we take the view of QNMs as modes that are smooth at
$\mathcal{H}^{+}$ and through a conformal compactification of
$\mathcal{J}^{+}$. The EFT correction does not affect the spacetime in these
regions, thus we do not expect it to change the BCs for QNMs. Because this
term, simply overloads notation, we will omit it in the remainder of the
discussion. Inserting this in equation (198b), we get:

$\begin{cases}\widetilde{\mathcal{O}}_{1,\omega}\widetilde{\Phi}_{y,x}^{(0)}=\widetilde{\lambda}_{(e)}\left(-\widetilde{\mathcal{O}}_{1,\omega}\widetilde{\Phi}_{y,x}^{(1)}+\frac{\alpha_{-1}}{\alpha_{1}}\widetilde{\mathcal{D}}_{\omega}\widetilde{\Psi}_{y,x}\right)+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\\\\[5.0pt]
\widetilde{\Phi}_{y,x}^{(0)}=\pm\,(-1)^{m}\widetilde{\mathcal{T}}_{0,\omega}\,\mathbb{P}\widetilde{\Psi}_{y,x}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{1})\end{cases}$
(201)

where $\,\widetilde{\,}\,\,$ identifies $\alpha_{\pm 1}$, accordingly, was
commuted past the corresponding operator, and we used
$\mathbb{P}\mathcal{\alpha}_{-1}=\mathcal{\alpha}_{1}$. Finally, we expand
$\omega$ in powers of $\widetilde{\lambda}_{(e)}$:

$\omega=\omega_{0}+\widetilde{\lambda}_{(e)}\,\omega_{1}+\mathcal{O}(\widetilde{\lambda}_{(e)}^{2})\,,$
(202)

where $\omega_{0}$ is the background QNM frequency. We substitute this in
(201) and equate each order of $\widetilde{\lambda}_{(e)}$ to $0$. The
$0^{th}$ order is satisfied by definition, and the $1^{st}$ order yields:

$L^{(0)}\,\widetilde{\Phi}_{y,x}^{(1)}=\left(\mp(-1)^{m}\omega_{1}\,\delta
L\,+L^{(1)}\right)\widetilde{\Psi}_{y,x}\,,$ (203)

where:

$\displaystyle L^{(0)}$
$\displaystyle=\widetilde{\mathcal{O}}_{1,\omega_{0}}\,,$ (204) $\displaystyle
L^{(1)}$
$\displaystyle=\frac{\alpha_{-1}}{\alpha_{1}}\widetilde{\mathcal{D}}_{\omega_{0}}\,,$
$\displaystyle\delta L$
$\displaystyle=\partial_{\omega_{0}}\widetilde{\mathcal{O}}_{1,\omega_{0}}\widetilde{\mathcal{T}}_{0}\,\mathbb{P}\,.$

Mapping $\omega_{1}\rightarrow-\omega_{1}$ takes the $+$ into the $-$
polarization. Thus, for each background frequency $\omega_{0}$, the two
possible corrections have opposing sign. This is in agreement with what was
postulated in equation (120) for $\lambda_{(o)}=0$. Consequently, we find
agreement with the results in the Eikonal limit, (section IV.2.1) and the
results derived for the Schwarzschild limit (section IV.2.2).
$\widetilde{\Psi}_{y,x}$ separates as the product of a radial and an angular
function. Thus, we can use the corresponding EOM to iteratively reduce the
order of the derivatives of $\widetilde{\Psi}_{y,x}$ in equation (203) into at
most first order. The resulting equation is in the form of equation (241b),
and can be solved with the pseudo-spectral methods outlined the in appendix A.

###  IV.3 Mixed parity corrections

In this section, we will generalize the discussion of IV.2 to mixed parity EFT
corrections, i.e. modes where both $\lambda_{(e)}$ and $\lambda_{(o)}$ are
non-zero. In section IV.2, the main non-trivial step was to find a basis such
that the EOM for the two possible electromagnetic polarizations decouple. In
the Schwarzschild limit, we found that this could be done by taking the real