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

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.

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 twoparameter Kerr family of metrics.
The mass and angular momentum (M, a) of the BH are sufficient to specify these highly symmetric solutions of the Einstein equations.In the limit a → 0, we obtain the only static, spherically symmetric solution of the equations of motions (EOMs), 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 (GWs) 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 (ODEs) [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 quasi normal modes (QNMs) 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 B 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].
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 ≲ 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 ≳ 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 (PDEs).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 Φ.We require that Φ is small and that, in the Φ → 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 Φ.Additionally, for simplicity, we restrict to terms that are parity invariant. 1fter suitable field redefinitions, the leading order EFT can be written as follows [36]: 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 Φ, we derive the EOM: Different choices of f (Φ) have different consequences.If f (Φ) is constant, the right hand side (RHS) vanishes, and we recover the massless Klein-Gordon (KG) equation.Conversely, if f (Φ) ∼ Φ, 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 (Φ) ∼ Φ 2 .For this theory, Φ = 0 solves the scalar EOM, and thus, GR BHs are in the solution space.Take f (Φ) = λΦ 2 /2 for some O(1) dimensionless coupling parameter λ, the EOM for the scalar field becomes: This equation is analogous to the massive KG equation − m 2 Φ = 0.In a flat background, when m 2 > 0, Φ 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 λ < 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 λ > 0, BHs with sufficiently small mass (M ) can be unstable.Indeed, the authors show that for certain mass ranges obeying λ L M 2 ≳ 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, G can be negative, reversing the sign of the effective mass.In this case, the Tachyonic instability emerges for λ < 0. Again, this only occurs for large enough |λ|(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 λ.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 To generate these plots, we computed the QNMs for BHs with a ≲ 0.999M in most cases, or a ≲ 0.999975M for ℓ = m ̸ = 0 modes.We found that, in general, the imaginary part of the background frequency increases with ℓ.Thus, the lowest lying QNMs should have ℓ → ∞.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 ℓ = m ̸ = 0 modes M ω0 → ℓ/2, as predicted in [38].This feature extends to the EFT correction: ω1 → 0 as a → M .
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: For the reader's convenience, we have summarized the main results in figure 1.We will discuss the results in detail in section III C. 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: FIG. 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 ℓ.Contrary to the scalar case, for most of the parameter space, Im(ω0) decreases with ℓ, making the mode with ℓ = 1 the lowest lying mode.This tendency is broken at a ≈ 0.96M , where the ℓ = 2 modes become the lowest lying QNM.As in the scalar case, for ℓ = m modes, in the limit a → M , we have M ω0 → ℓ/2.This is in agreement with the prediction in [38].Furthermore, as seen in the scalar case, the corresponding correction obeys M ω1 → 0.
Here, L is the EFT length scale, F represents the electromagnetic tensor, C denotes the Weyl tensor, ⋆ denotes the hodge dual operation (see section II F), and λ (e) /λ (o) are dimensionless coupling constants.The terms proportional to λ (e) and λ (o) are, respectively, even / odd under parity transformations, see section II E 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 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.
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 λ 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. where: 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 ω The corrected QNM frequencies are not centred around ω 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 ≳ 0.75 that would be difficult to capture via a perturbative approach.We will discuss these results in detail in section IV D 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 λ (e) is known [45]: (9) Here, M p is the Planck mass, m e is the mass of an electron, and M ⊙ 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 λ (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, • • • for abstract index notation, Greek indices µ, ν, σ, • • • for quantities in a coordinate basis, and parenthesized Greek indices (α), (β), (γ), • • • for quantities in a null tetrad basis.

A. 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, ϕ), that differ from standard BL coordinates by the definition x = cos θ.The line element reads: where: We will restrict our study to BHs such that M ∈ R + .In the limit a → 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 , ∆ 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. H + and H − are the BH and white hole horizons, respectively, whereas J ± 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 (CH).The surface gravity κ of H + is: where: This parameter describes the angular velocity of H + as measured by an observer at infinity.In the extremal limit, r − → r + , and κ → 0. In this limit, we no longer have a Cauchy horizon, and H becomes a doubly degenerate Killing horizon.This behavior makes calculations in the extremal limit increasingly difficult.
Finally, we define the dimensionless angular momentum: If we replace a by J in the metric, a change M → αM is equivalent to the coordinate transformation r → r/α , t → t/α.Consequently, all dimensionless background quantities must be independent of M .J is the only physically relevant parameter.
Penrose diagram for the exterior of a Kerr BH.

B. 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 (α) .To define this structure, we begin by selecting two null curve congruences, γ l and γ n , with tangents l a and n a satisfying l•n = −1.Subsequently, we choose two spacelike vector fields X a and Y a orthogonal to both l a and n a , such that: These are then used to define complex vector fields m a and ma : Throughout this paper, the bar (¯) operation will denote complex conjugation.As a direct consequence of equation (15), we find m 2 = m2 = 0 and m • m = 1.Finally, the null tetrad will be e a (α) := (l a , n a , m a , ma ) (α) .The normalization conditions reduce to: with: After fixing the directions of γ l and γ n , there are still multiple consistent choices of e (α) , 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 γ l and γ n , such that: 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 θ, which, again, need not be a constant: The two transformations can be encoded in a single scalar field λ = √ re iθ/2 : In general, if under (21), η transforms as: we say that η is GHP-covariant, with type (or weight) p, q.We deduce that the tetrad components have type: To preserve GHP covariance, we cannot add two quantities with different weights.Conversely, multiplying a {p, q} scalar by a {p ′ , q ′ } scalar is allowed, yielding a {p + p ′ , q + q ′ } scalar.To further simplify the notation, we introduce the prime ( ′ ) operation.This effectively interchanges instances of l a and n a , and swaps instances of m a and ma : If η is of type {p, q}, then η ′ will be of type {−p, −q}.
Conversely, η will be of type {q, p}.Both operations commute, and (η ′ ) ′ = η = η.It is also useful to define the spin and boost weights of η as s = 1 2 (p − q) and p = 1 2 (p + q), respectively.Under complex conjugation, the spin weight's sign is flipped, whereas the boost weight remains unchanged.Under a ( ′ ) 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 groups2 : and together with the corresponding primed versions.One can show that the scalars in (25) are GHP-covariant, with types: Conversely, because λ is spacetime varying, the scalars in (25) do not transform according to (22), gaining an additional affine contribution: A similar problem arises when taking the covariant derivative of GHP-covariant quantities.An affine contribution will stem from the derivatives of λ and λ.Consider l a ∇ a A (l) := l a ∇ a l b A b .Under (21) we obtain: Comparing equations ( 28) and ( 29), we deduce that for the quantity (l a ∇ a − ε − ε) 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: where the numbers on the right denote the GHP weight.þ, þ ′ respectively raise and lower the boost-weight, leaving the spin-weight unchanged.Conversely, ð and ð ′ 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: This operator obeys a Leibniz rule.For α, β GHPcovariant quantities: 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: these are usually known as the NP scalars.Similarly, the 6 components of the Maxwell tensor, can be encoded in 3 complex scalars: with 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]: 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: Thus, the only non vanishing GHP scalars are ψ 2 , ρ, ρ ′ , τ and τ ′ .The Kinnersley tetrad is an example that satisfies this condition [43,57]: C. 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]: where and we defined Θ a in equation (34).When acting on a {0, 0} qunatity, Θ a reduces to ∇ a , thus T 0 is the standard d'Alembert operator .The Teukolsky master equation is a generalization of this that acts on ψ (s) , a spin-s field with weight {2s, 0} (see table I).The master Teukolsky equation reads [58]: Using the Kinnersley tetrad (equation ( 40)), this equation can be written: Now, to separate the resulting PDE, we follow the approach in [3].We make the Ansatz : where and m ∈ Z, ensuring periodic boundary conditions (BCs) are met for ϕ.Substituting this in the Teukolsky equation, we obtain: where O s,ω is an operator that only depends on (r, x).
Finally, taking Ψ s = R s (r)Θ s (x), we can separate equation (47) as the sum of a radial and an angular ODE operator [3,16,59]: where: with: and s A ℓm is a separation constant.
In the ψ4 (1)   TABLE I. Physical meaning of the fields that solve equation (43).
Schwarzschild limit, regular solutions of equation ( 49) are the spin-weighted spherical harmonics, with s A ℓm → l(l + 1) − s(s + 1), see [60].When we turn on rotation, the equation depends on ω, 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 = ±1 for (49) and r = r ± 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].

D. 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 H + and J + , see [62][63][64].Physically, this is the statement that QNMs are waves that are ingoing through H + and outgoing at 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]: where: with: Similarly, an asymptotic expansion near r → ∞ yields a decomposition into an ingoing and an outgoing mode: where: QNMs are defined as solutions of the Teukolsky equation with A(ω) = D(ω) = 0.This defines a boundary value problem (BVP) with eigenfunction R s and eigenvalue ω.
Numerically, there is a more practical way of setting up the BVP.Define the compact radial coordinate: At y = 0 (r = r + ), and y = 1 (r = +∞), the BCs corresponding to ingoing and outgoing modes and the multiplicative factor that transforms one into the other are singular.Define: where α y,s encodes the singular part of the appropriate BCs, [16]: Substituting this in equation ( 50), and requiring smoothness of R s (y) in y ∈ [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 ∈ [0, 1], it must describe a QNM.
In modified BL coordinates, the Kerr metric has coordinate singularities at both S 2 poles (x = ±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: with: QNM solutions will have Θ s (x) to be a smooth function 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 α 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 Teukolsky modes, the joint singular behaviour is: Substituting this in equation ( 47), we get: QNMs are solution of the full PDE that are smooth for

E. Parity transformations
Denote P the parity transformation around r = 0.In modified BL coordinates, P takes x → −x and ϕ → π+ϕ.The action ( 6) is invariant under parity, and so is the Kerr metric.For the Kinnersley tetrad,see (40), we can prove: 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 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 → m and m → m.This motivates the definition of the conjugate-parity-transform operator P. For some scalar η, we have: As an example, consider the action in the Teukolsky field ψ ( s).Using equation ( 45), we have that: where the (−1) m contribution comes from taking ϕ → π + ϕ.
The transformation rule in equation ( 65) is tensorial.Thus, if η is constructed by taking covariant derivatives and tensor products of tetrad quantities, we have: where s is the number of m a and ma instances in the definition of η.As an example, take ρ and τ .Using equation ( 25), we have: By definition, η is a GHP covariant quantity.Thus, s is simply its spin weight: s = (p − q)/2.Now, consider 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: where s is the spin-weight of G.Because Kerr is invariant under parity, this rule also extends to the case where G contains pure gravitational scalars, e.g.ψ 2 .
In particular, the Teukolsky operator O s has spinweight 0, thus: where [•, •] denotes the commutator operator.Standard results tell us that if ψ (s) is in the kernel of O s , then so is Pψ (s) .Using equation ( 45), we have that: This proves that if (m, ω) is in the QNM spectrum of O s , then so is (−m, −ω).The parity invariance of the Kerr metric leads to a degeneracy between these two families of QNMs.
The action of 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: Thus, to preserve (65), we must have: The absolute value and phase of λ must be parity even and odd respectively.Take η to be a GHP scalar with boost weight r and spin weight s.Restricting to GHP transformations that respect (74), we get: The action P is akin to complex conjugation.It takes a {p, q} GHP scalar to a {q, p} scalar.Consequently, P preserves the GHP weight.This fact is particularly useful, as it allows us to take linear combinations of η and Pη without breaking GHP covariance.

F. 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, ⋆, is a linear map from the space p forms to the space of d − p forms.We will follow the convention: where ε is the Levi-Civita tensor.We have For a Lorentzian metric, we have: Finally, we can define the differential and co-differential operators: 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 C: Now, taking the left and right dual at the same time, we get: (⋆C⋆) In the second step we used the standard result for the contraction of two instances of ε, and in the the third step, we used the fact that C is traceless, i.e.C c acb = 0. Using ⋆ 2 = −1 (for a 2 form in 4 dimensions), we get: 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 ⋆, we decompose the electromagnetic and Weyl tensor into their self-dual and antiself-dual parts: As intended, F ± and C ± are ⋆ eigenstates: When projected to the Kinnersley tetrad ( 40), the Levi-Civita tensor takes a simple form: where E is fully antisymmetric and E 1234 .Thus, we get Dualized quantities, carry an extra factor of i, and thus, complex conjugation will add an extra − contribution.Concretely, consider the dualized GB term: 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 E. Thus, under the conjugate-paritytransform, we have: As a rule of thumb, quantities that are dualized an odd amount of times are parity odd.

G. 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 ϕ 0 , and then explain how to generate F ab from it.
The vacuum Maxwell equations take the form: 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 F, to combine the two equations into a single complex equation: 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]: The corresponding GHP weights are: If we increase the spin-weight of M l and the boostweight 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 ζ introduced by Teukolsky in [3]: Using the GHP relations in [53], we can obtain the Teukolsky equation for s = 1 (see equation ( 44)): F is a closed form, thus, locally, there is a form A, such that F = dA.This is known as the vector potential.We can express the components of F directly in terms of the A. Using the definition of ϕ i (equation ( 37)), and treating A a as a {0, 0} GHP covariant vector field, we have: Thus, in terms of A, the Teukolsky equation is: On the other hand, working directly in an abstract basis, we can obtain M a in terms of A a .We have: Acting with ζ a on the left, we must recover the Teukolsky equation: Thus, we have the operator identity: The key insight comes from taking the transpose4 of this identity.We have that: Now, we can prove that E T = E, thus, for a mode ψ H such that: we have: We proved that Âa = ζ a T ψ 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, ψ H is a complex scalar, thus, Â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).E is a real linear operator, thus if Â is in its kernel, then so must be Â.Hence, for any mode ψ H that solves (99), we can construct a real solution to Maxwell's equation given by: In the literature, ψ H is usually known as the Hertz potential.
From equation ( 43), we can deduce that O s T = O −s .Thus, ψ H is a solution to the master Teukolsky equation with s = −1.In fact, this is a re-parametrization of the EOM for ϕ 2 [58]: Thus, the Hertz potential is proportional to the ϕ 2 component of a solution to the Maxwell equations, and we can deduce ψ H has GHP weight {−2, 0}.We can prove that any QNM of ϕ 2 will generate a QNM for A a with the same frequency, and consequently, all QNMs of ϕ 2 must be QNMs of ϕ 0 .The Hertz potential derivation can be reversed, and we can generate a vector potential A a from any ϕ 0 .Thus, all ϕ 0 QNMs must also be QNMs of ϕ 2 .Hence, ϕ 0 and ϕ 2 must be isospectral, as seen in [43].
Notwithstanding, the ϕ 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 τ i defined in (II G).Applying the GHP identities in [53] we get [65]: Because we required A to be a real vector, φi can be obtained by simply taking the complex conjugate of these expressions.If we worked with complex A, ϕ i and φi would be independent variables.Conveniently, ϕ i depends exclusively on ψH . 5Thus, we can encode equation (103), in operators T i such that:

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 B, we will briefly outline how to adapt our methods to consider parity odd corrections.

A. 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 C, they reduce to: where O 0 is the s = 0 Teukolsky operator and G is the GB term, defined in (2).Because Kerr is a vacuum, type-D spacetime, we get: Following the discussion in section II E, it is clear that the conjugate-parity-transform operator, P, commutes with all terms in (105): thus, if Φ is a solution to the EOM, then so is PΦ.
As argued in the end of section II E, Assuming Φ is proportional to S(t, ϕ) = exp(−iωt + imϕ), we get that if (m, ω) define a QNM, then so does (−m, −ω).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, Thus, the EOM are, in general, not separable.Notwithstanding, we can still follow the procedure in section II D to define the QNMs.First, as seen in ( 45), we factor out the (t, ϕ) dependence, by making the Ansatz : We get: where O 0,ω was defined in equation (47).Now, we must enforce BCs corresponding to outgoing behavior at J + and ingoing at H + .To do so, we will replace r with the compactified radial coordinate y = 1 − r/r + .In the case λ = 0, as argued in section II D, we can set the BCs by factoring out α 0 (r, x) (see equation (62)) and requiring the remainder to be smooth for (y, x) G is regular at r + and at x = ±1.Furthermore, it is subleading with respect to O 0,ω at r → ∞.This implies that the leading order of a Frobenius expansion at the domain boundaries must be independent of λ.Thus, QNMs of the corrected equation obey the same BCs as the QNMs in the uncorrected case.Define: Replacing this in equation ( 110), we get: QNMs are solutions of this equation, that are smooth for (y, In the Schwarzschild limit, the RHS reduces to 48L 2 M 2 λ/(r + (1 − y)) 6 Ψ(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 ℓm .The radial part is an eigenvalue problem for ω, with (L/M ) 2 λ and ℓ 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 ̸ = 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 λ = O(1), the relevant perturbation parameter is: Our results are only valid to leading order in λ.
We can solve the EOM using standard quantum mechanical perturbation theory.We make the following Ansatz : Grouping together same powers of λ, and neglecting higher order terms the EOM reduce to: 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 ( Ψ (1) , ω 1 ).
Apart from ( Ψ (1) , ω 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.

B. 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 G with: As argued at the end of section II F, P anticommutes with G: thus, the action of P in the EOM will be: P flips the sign of λ.Thus, again, we can relate QNMs whose frequency has positive real part with the ones with negative real part.If (m, ω(λ)) defines a QNM frequency, then so does (−m, −ω(−λ)).Bearing this in mind, and following the approach in section III A, we could obtain the parity odd EFT corrections to scalar QNMs numerically.

C. 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 ω 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 ω 1 increases exponentially.
The main results are presented in figure 1 and table II.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 ℓ = m = |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 ℓ → ∞.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 ℓ 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 ℓ.An in-depth discussion of this lies beyond the scope of this paper.FIG. 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 ∈ (0, 0.999975), by solving equation ( 116).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.ℓ = m ̸ = 0 modes vary very quickly near extremality, with M ω → ℓ/2 as a → M .If we truncated the spin at a/M ∼ 0.99, we would lose roughly the last third of the parametric curves.FIG. 6.In this image we explore the scalar QNM frequencies of ℓ = m ̸ = 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 → M , M ω0 → ℓ/2, agreeing with the prediction in [38].Similarly, in this limit, M ω1 converges to 0. Notice that the near-extremal behaviour of these modes has some structure.In particular, M Re(ω ℓ=1 0 ) crosses 0 at a/M ≈ 0.9993, and M Im(ω ℓ=2
In the Schwarzschild limit, the spacetime is spherically symmetric, and the EOM can be solved nonperturbatively in λ, using the direct methods in section A 2. We obtained the QNM frequency ω for several values of λ.Taking λ to be arbitrarily small, we disentangled ω into a background component ω 0 and the EFT correction ω 1 .The results can be seen in the plot markers of figure 5.For a ̸ = 0, the EOM are no longer separable, thus we proceed perturbatively, as outlined in section III A. The results are shown in the plot lines of figures 1 and 5, and in table II.
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 ω 1 have maximal absolute value for ℓ = m modes.For ℓ ̸ = 0, the correction is largest at a ∼ 0.9M , quickly approaching 0 afterwards.For a ≳ 0.9M , the correction to modes with ℓ ̸ = m > 0, or ℓ = m = 0 increases quickly, before approaching finite values at a → M .This can more easily be seen in a parametric form in figure 5. ℓ = m ̸ = 0 modes, also known as FIG. 7. Here, we show M Im(ω0) for scalar modes in the Schwarzschild limit as a function of ℓ.The dark line represents the Eikonal limit prediction.Contrary to gravitational QNMs, Im(ω1) increases with increasing ℓ, approaching the Eikonal limit frequency from below.This indicates that the longest-lived modes will have ℓ → ∞.
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 ≲ 0.999975.As predicted in [38], we find that ω 0 approaches ℓ/2.Furthermore, for ℓ = m = 2, the convergence rate of M ω 0 − 1 seems to be consistent with 1 − a/M , as predicted in [38].However, for ℓ = m = 1, M ω 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 ∼ 0.999975M .On the other hand, we found that, for slowly damped modes, M ω 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 ω 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 cer-tain 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 ℓ 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 B, we will set λ (o) = 0 in the action (6), and obtain the corrections to the two families under parity even EFT corrections.Then, in section IV C, we will keep λ (o) ̸ = 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: where: and: The factor of (−1) ℓ+m+1 was introduced to ensure that, in all cases studied, the Schwarzschild limit of ω 1 has positive real and imaginary parts. 6Because ω + 1 = −ω − 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 D, we will use the pseudospectral methods in the appendix A to solve the EOM, obtaining the corrections to the QNM frequencies in all cases.

A. 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: We want to solve this equation directly in terms of the components of F ab , thus, we must also enforce the electromagnetic Bianchi identity, dF = 0, guaranteeing F is a closed form.Following the approach in section II G, the two conditions can be encoded by the real and imaginary part of a single complex equation.We get: where, F − was defined in equation ( 82), and we define ":" as a double contraction.For some 4 tensor A ab cd and 2 tensor B ab we have: Assuming 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 λ (e) and λ (o) into a single complex coupling constant: Using the identities in section II F, we get: By definition, F + = F − , thus the EOM depend on both, F − and its complex conjugate.Consequently, the equation is not holomorphic in F ± , and we must independently solve its real and imaginary parts.
To make progress, define F ± , such that: Note that F ± obeys the same properties as F ± .We have: Thus, we can define: as a 2 form that reduces to F in the L → 0 limit.Working perturbatively in L/M , in terms of F ± , the EOM are: To project this into a null tetrad basis, we must first define NP scalars to encode F. Following the convention in equation (37), we get: (132) In the limit L → 0, Φ i → ϕ i .The Weyl tensor also obeys a Bianchi identity ∇ [e C ab]cd = 0.In a vacuum, this implies that ∇ a C abcd = 0. Putting the two identities together, we can prove: Using this, we can prove that in the null tetrad basis, (131) reduces to: where: and: .
As expected from the identities n a = (l a ) ′ and ma = (m a ) ′ , the 2 nd and 4 th line can be obtained from the 1 st and 3 rd by using the prime operation.
While M (α) depends exclusively on Φ i , J (α) , is a function of Φi .Thus, if we try to factor out exp(−iωt) from each Φ i , we will get a term proportional to exp(−iωt) on the LHS and a term proportional to exp(iω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 Φ i and Φ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 B, 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 C.

B. Parity even corrections
In this section, we set λ (o) = 0 and consider only parity even EFT corrections.In section IV B 1, Using the effective metrics derived in [22], we will obtain an analytic expression for the QNMs in the large ℓ limit.Then, in section IV B 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 B 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.

Eikonal limit
The Eikonal approximation, also known as the geometric optics limit, is valid when the wavelength of light, dubbed ε, is much smaller than the characteristic length scale of the background spacetime.For a Kerr background, this is the statement that α ≪ 1, with: 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 ℓ limit of QNMs to derive analytic expressions for the frequencies in this limit.
We start with a plane wave Ansatz : We assume that the characteristic length scale of all variables is M , with the factor of 1/α in the exponential ensuring that the signal has a short wavelength.Plugging this into equation ( 123), and setting λ (o) we get: where: K a := ∇ a Θ 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 C ∼ 1/M 2 , thud, the size of the correction in (140) is: Now, the EFT approximation is only valid, if L is much smaller than any other length scale in the theory.In particular, L/ε ≪ 1, or else, higher order derivative corrections become relevant.Thus: In the Eikonal approximation, we neglected all terms that are O(ε/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 ℓ limit.Roughly, we have: The angular part of a QNM with azimuthal number ℓ has ℓ nodes, thus, the length-scale of this mode is at least M/ℓ.Hence, the validity of the EFT approximation requires L/M ≪ 1/ℓ.We expect the correction to the QNM frequency to be O(L/M ) 2 smaller than the background, thus it should only be present at O(1/ℓ 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 ℓ QNMs, ignoring the requirement L ≪ ε.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 A a is in the same direction as K a , however, this is a pure gauge solution (A a ∼ ∇ a e iΘ ), 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: where 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]: where: Here, g ab is the background Kerr metric, and: For L ≪ M , g (±) ab are non-degenerate symmetric 2tensors, 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 D 2, we will work perturbatively in L/M .At leading order, we have: Substituting this in equation ( 146), we get: where: For a Kerr spacetime, the large ℓ limit of QNMs is fully specified by circular photon orbits [13][14][15].In fact, ℓ = ±m QNMs are related to the two sets of unstable equatorial null geodesics.The following relation holds: (151) Here, Ω = dϕ dt represents the angular frequency of the orbit, and γ 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.ℓ = m modes relate to null-geodesics that co-rotate with the BH, while ℓ = −m modes are specified by counter-rotating photons.
Perturbatively, the two effective metrics are related by reversing the sign of λ.Thus, in the Eikonal limit, the corresponding QNMs must be related by reversing the sign of λ.Thus, following the convention in equation (120), we can conclude that: Null geodesics extremize the Lagrangian, where χ is some affine parameter.By definition L = 0, furthermore due to the symmetries of g, we have that: are conserved quantities.Inverting this relation, we can express (dt/dχ, dϕ/dχ) in terms of (E, L).Because we are looking for equatorial geodesics, x = dx/dχ = 0. Substituting all of this in L = 0, and solving with respect to dr/dχ, we get: where: and: is the impact parameter of the null geodesic.Taking the derivative of (155) with respect to χ, we see that circular geodesics happen if and only if: To solve this, we will make the Ansatz : At O λ (±) 0 , the equations reduce to: FIG. 8. Here, we show the prediction for the large ℓ 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 ℓ = −m modes, with a/M ∈ [0, 1], running from the middle to the left.Conversely, the right side illustrates the ℓ = m modes, with a/M ∈ [0, 1], from the middle to the right.We arrange the plots in this manner to underscore the equivalence between the m → −m operation and a → −a.In the near extremal limit, for ℓ = m modes, we find that M ω ℓ=m 0 → ℓ/2, and M ω ℓ=m 1 → 0. Finally, note that in the Schwarzschild limit, Im(ω1) → 0, as can be deduced from equation (165) The cubic equation has two roots for r > r + , given by [14]: Notice that r + 0 ≤ 3M ≤ 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 ± is independent from the one in λ (±) , 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: We can now compute the relevant parameters for equation (151).The Lyapunov exponent is defined as: Expanding in powers of λ (±) , we get: with: Similarly, the orbital frequency is given by: 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 ℓ limit.
The results are summarized in figure 8.

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 nonperturbatively 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 ≪ 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: Furthermore, all scalars depend exclusively on the radial coordinate r.Projecting equation (127) into a tetrad basis, and setting λ o = 0, we get: To simplify these equations, our strategy is to express them as a system of real equations of real variables.ϕ 1 has spin weight 0 (see section II B), 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 ϕ 1 , without breaking GHP covariance.This would not be possible for ϕ 0 and ϕ 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 ð and ð ′ respectively.Define the polarized Maxwell GHP scalars as: 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 spinweight 1 and −1 respectively.Thus, after applying ð ′ and ð 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: and: In the Schwarzschild limit, þ and þ ′ reduce to (t, r) derivative operators, while ð and ð ′ act exclusively on the angle variables.Thus, the only angular dependence in equations ( 170) and ( 171), is encoded in the (ð ð ′ ) operator.There is a known relationship between (ð ð ′ ) and the spherical harmonics Y ℓm [55].We have: Given þY ℓm = þ ′ Y ℓm = 0, a decomposition of the scalars in spherical harmonics removes all angular de-pendence.Similarly, we can factor out time dependence by multiplying out exp(−iωt).7 P ± 0 and P ± 2 have boost weight 1 and −1 respectively.Thus, to work with {0, 0} scalars, we will explicitly factor out ρ and ρ ′ respectively.The final Ansatz will be: 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 (r) and dR ± 1 (r)/dr from the 3 rd and 4 th equations, reducing the problem to a system of two ODEs with two variables.Finally, we define u ± (r) and v ± (r), such that: with where we defined: The EOM reduce to: where: These equations can be combined into a single decoupled second-order ODE for v ± .We get: After suitable variable changes, we can prove that this equation agrees with the effective potentials derived in [51] (equations 16 and 17). 8In the λ (e) → 0 limit, Ξ → 1, ζ → 1, and the EOM reduce to the wave equation for electromagnetic fields in a Schwarzschild BH background, see [67].
As outlined in section II D, to find the QNM spectrum of (179), we must set BCs that are ingoing at r → r + , and outgoing at r → ∞.Ξ and ζ are smooth at r = r + , and asymptote to 1 + O(1/r 3 ) at r → ∞.Thus, at leading order, a Frobenius expansion near the boundaries is independent of λ (e) .Hence, to set the appropriate BCs, we can proceed exactly as in section II D. Define v ± (r): with: 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 ∈ [0, 1].Treating λ (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 λ (e) , and can, in theory, be solved for λ (e) ∈ (−12, 6). 9However, the EFT approximation entails λ (e) ≪ 1, thus, we should disregard the remaining of the parameter space.Expanding equation (179) in powers of λ (e) , we see that at leading order, the EOM for the ± polarizations can be mapped into each other by changing λ (e) → − λ (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). 10

Generalization to Kerr
In this section, we will solve equation ( 134), considering only parity even EFT corrections (i.e.λ (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 λ (e) , defined in (176).
Take a generic differential equation: where Q (0) , G and Q (1) are arbitrary.At O( λ 0 (e) ), we have ). Substituting this on the O( λ 1 (e) ) part of the equation, we get: Thus, any identity that is satisfied at O( λ 0 (e) ) can be used to simplify the O( λ 1 (e) ) terms.Following this idea, we equate M (α) = 0 and solve with respect to the GHP derivatives of ϕ 1 .The result can be used to eliminate all derivatives of φ1 from J (α) .Then, we apply the Teukolsky ζ a operator, defined in equation ( 91), to both sides of equation ( 134).Using the GHP identities in [53] we get: where the summation in i is implicit, J i are linear operators: and O s is the master Teukolsky operator defined in equation (44).
From section II G, we know that there is a Hertz potential ψ H such that: where T i was defined in equation ( 104).Substituting this in equation ( 184), we get: with: Equation (187) mixes modes proportional to S with modes proportional to S. If we set Φ 0 and ψ H proportional to S = exp(−iωt + imϕ), 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 S on the RHS.The reverse problem occurs if we choose ψ H to be proportional to S instead.In section II E, we studied the properties of the conjugate-parity-transform operator, 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 O s .Consequently, we found that if ψ (s) = S(t, ϕ)Ψ s (r, x) is a QNM of O s then so is Pψ (s) = (−1) m S(t, ϕ)Ψ s (r, −x).Taking linear combinations of these states, we can express the kernel of O s in a basis of P eigenstates.We restricted to parity invariant EFT corrections, thus, we should be able to express solutions to equation (187) as eigenstates of P. Take: By definition, PΦ ± 0 = ±Φ ± 0 and Pψ ± H = ±ψ ± H , thus we associate these with the polar and axial polarizations of light respectively.11Because, O s , D and T 0 have even spin-weight, they all commute with P: Thus, applying P to both sides of equation ( 187), we get: where P is the parity transform operator, (see section II E).Writing equations ( 187) and (191) in the P eigenstates basis, we get: Because [P, O −1 ] = 0, we have O −1 ψ ± H = 0.However, in a Kerr background, QNMs of the Teukolsky equation are not, in general, eigenstates of P, consequently ψ ± H cannot be a QNM solution.The same can be said about the O( λ 0 (e) ) component of Φ ± 0 .Nevertheless, we can choose ψ ± H and Φ ± 0 to be a linear combination of modes proportional to S and S. Extending this to first order in λ (e) , we take: with Ψ and Φ proportional to S:12 PΨ and PΦ will be proportional to S, thus substituting this in equation ( 192) we get: where: S and S are linearly independent, thus, their pre-factors must vanish.We get: Equation (198b) relates the background of Φ r,x with the background of Ψ 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 λ (e) → 0, Φ must be a QNM of the s = 1 Teukolsky equation.We substitute this in equation (198a), and work perturbatively in λ (e) to find the QNM spectrum of the full theory.Concretely, we will solve (198a) in two steps.First, treating ω as a parameter independent of λ (e) , we set the appropriate QNM BCs.As outlined in section II D, this is done by factoring out the appropriate corresponding singular behaviour.Then, we expand ω in powers of λ (e) , and solve equation (198) 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 Ψ r,x is well understood.Working with the compactified radial coordinate, we have: where α −1 , defined in (62), is singular at the domain border, and Ψ y,x is smooth for (y, x) Working perturbatively in λ (e) , Φ r,x takes the form: Φ r,x (r(y), x) = α 1 (y, x) + λ (e) δα(y, x) Φ (0) y,x (y, x) + λ (e) Φ (1)  y,x (y, x) By definition, Φ . For Φ to be a QNM, we require Φ (1) y,x to be smooth in the same interval.δα 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 δα = 0. 13 Because this term, simply overloads notation, we will omit it in the remainder of the discussion.Inserting this in equation (198), we get: where identifies α ±1 , accordingly, was commuted past the corresponding operator, and we used Pα −1 = α 1 .
13 This is expected if we take the view of QNMs as modes that are smooth at H + and through a conformal compactification of J + .The EFT correction does not affect the spacetime in these regions, thus we do not expect it to change the BCs for QNMs.
Finally, we expand ω in powers of λ (e) : where ω 0 is the background QNM frequency.We substitute this in (201) and equate each order of λ (e) to 0. The 0 th order is satisfied by definition, and the 1 st order yields: where: Mapping ω 1 → −ω 1 takes the + into the − polarization.Thus, for each background frequency ω 0 , the two possible corrections have opposing sign.This is in agreement with what was postulated in equation (120) for λ (o) = 0.
Consequently, we find agreement with the results in the Eikonal limit, (section IV B 1) and the results derived for the Schwarzschild limit (section IV B 2). Ψ 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 Ψ y,x in equation ( 203) into at most first order.The resulting equation is in the form of equation (A1), and can be solved with the pseudo-spectral methods outlined the in appendix A.

C. Mixed parity corrections
In this section, we will generalize the discussion of IV B to mixed parity EFT corrections, i.e. modes where both λ (e) and λ (o) are non-zero.In section IV B, 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 and imaginary parts of suitably transformed EOM.These components are invariant under complex conjugation, thus, by taking the real and imaginary part of the EOM, we are effectively projecting them into a basis of complex conjugation eigenstates.To generalize the result to Kerr BHs, we followed a similar approach, projecting the EOM into a basis of eigenstates of the parity complex conjugation operator P. In both cases, we expected the approach to work because of the parity invariance of the equations of motion. 14 In this section, we are no longer dealing with parity invariant EFT corrections, thus the approach should break down.Notwithstanding, it may be possible to choose a different basis that shares the same properties as the ones in section IV B, where the equations decouple.
Define A δ , with δ ∈ [−π, π), as an operator that takes α ∈ C and projects perpendicularly onto the line in the complex plane defined by re iδ , r ∈ R. We get: As expected, A 0 α = Re(α) and A π 2 α = Im(α).Furthermore, note that the output is always a real number, and thus invariant under complex conjugation.Crucially, if we know the projection of α at two different angles, δ 1 , δ 2 , as long as they do not differ by an integer multiple 14 In the Schwarzschild limit, due to spherical symmetry, when we take the real and imaginary part of the EOM, we are also indirectly projecting into a basis of eigenstates of P. We did not make this more explicit because the EOM look slightly more complicated.
of π, we can always recover α and ᾱ: A δ , allows for a systematic way to find a basis to express (α, ᾱ) whose components are real valued.Note that A δ is periodic in δ.For any integer k, we have: Furthermore, flipping the sign of δ is equivalent to taking the complex conjugate of the argument: Finally, note that rotating α by some phase γ ∈ [−π, π), is the same as rotating the δ line in the opposite direction, i.e. δ → δ − γ: Now, lets employ this approach to a toy problem, similar to equation ( 131).Imagine we want to solve the following complex valued differential equation: Here, A, B are real valued differential operators, ψ is a complex valued function, and λ ∈ C has phase λ θ ∈ [−π, π): Because the equation depends on both ψ and ψ, it is not holomorphic in ψ.Thus, we need to explicitly solve the real and imaginary parts of the equation, and this is in general non trivial.However, there might be a lines in the complex plane where the equations decouple.Applying A δ to both sides, and using the fact that A and B commute with complex conjugation, we get: Now, we want to pick two distinct values of δ such that (212) depends on a single real function, i.e.A δ ψ ∼ A λ θ −δ .Using equation (207), we get two distinct classes of values for δ, labelled as δ ± : where k ± ∈ Z. Solving this, we get: Plugging this into (212) we get two decoupled equations: We can solve this in terms of A δ+ ψ and A δ− ψ 15 , and then, using equation ( 206), recover the value of ψ and ψ.Note that the second equation can be obtained from the first by mapping |λ| → −|λ|.
Working perturbatively in L/M , in section IV C 1, we will use this approach to correct the QNMs of Schwarzschild BHs.Then, in section IV C 2, we will generalize the approach to find the corrections to QNMs of Kerr BHs of any spin.To do so, we will define an analogue of A δ that depends on P instead of the complex conjugation operator.In both cases, we will show that QNM frequencies reduce to the results of section IV B with λ (e) replaced by |λ| = λ 2 (e) + λ 2 (o) , in agreement with equation (120).

Schwarzschild limit
In this section, we will derive the EOM for the corrected QNMs of a Schwarzschild BHs, when considering mixed parity EFT corrections.In section IV B 2, we computed the corrected EOM for λ (o) = 0, non perturbatively in L/M .It is much harder to proceed non-perturbatively for mixed parity corrections, as there isn't an obvious way to decouple the EOM.Nevertheless, validity of the EFT approximation requires L ≪ M .Thus, only the leading order correction is physically relevant and we can safely proceed perturbatively.
Our starting point is equation (134).Making use of the GHP identities that are specific to spherically symmetric backgrounds (equation (167)), the EOM simplify.We get: where: and: 15 All values of k + , k − yield the same EOM.
The EOM are non-holomorphic, as they depend on both Φ i , and Φi .To tackle this, we will project the equations along suitably chosen directions of the complex plane, as outlined in the beginning section IV C.However, the action of A δ breaks GHP covariance.When we take A δ Φ 0 , we are adding Φ 0 and Φ0 , GHP scalars with different spin weights.Thus, as was done in section IV B 2, we must reduce / increase the spin weight of Φ 0 /Φ 2 respectively, to work with spinless scalars.We define the following set of quantities with 0 spin weight: Similarly, we must work with EOM of vanishing spin.To set the spin weights of the m and m components of (216) to 0, we act with ð ′ and ð respectively.The resulting EOM will contain factors of ðð ′ P 1 and ðð ′ P1 .As argued in section IV B 2, there is a relationship between the ðð ′ operator and the spherical harmonics Y ℓm , see (172).Using the fact that in the Schwarzschild limit, þ, þ ′ , are exclusively (t, r) derivatives and ð, ð ′ act exclusively in (x, ϕ), we can generalize this relation.For η, a GHP scalar with spin weight 0, that can be decomposed as a linear combination of spherical harmonics with angular quantum number ℓ: we have: Complex conjugation preserves the angular quantum number, thus it is safe to assume that P 1 and P1 obey this condition.In terms of P i the 0 spin-weight EOM are:16 where the summation in i is implicit and: Note that the LHS depends exclusively on P i , whereas the RHS is a function of Pi .Furthermore, for all i, M i and J i are real operators, the only complex parameter in (222) is λ.Thus, the EOM have the same form as (210), and we can decouple them into two polarizations following the same procedure. 17We find two phases, δ ± , that decouple the EOM.Taking λ θ to be the phase of λ, we get: The decoupled variables are: And the EOM separate as: Now, to derive radial decoupled EOM, we follow an identical procedure to the one in section IV B 2. First, we factorize the angular and radial dependence: 18 The spherical harmonics factor out, as there are no angular derivatives in (226).For each polarization, we get a system of four first-order ODEs of three radial variables.This is over-determined, thus, as was done in IV B 2, we can find linear combinations of the first two equations that eliminate R ± 1 (r) and dR ± 1 (r)/dr from the 3 rd and 4 th equations.The remaining two equations can be decoupled by setting: where: (229) 17 Note that in (210) we work with a single scalar ψ, whereas here we have a vector ψ i .Notwithstanding, because all J i commute with complex conjugation, the same approach works. 18P ± i are real variables, but we chose a complex valued Ansatz.Because the ODE operators are real valued, we can fix this by taking the real value of the solution in the end.
The EOM reduce to: where: Expanding equation ( 179) to leading order in λ (e) , and replacing λ (e) with (L/M ) 2 |λ|, we recover (230).Thus, to leading order in L/M , we find agreement between the QNMs obtained directly for parity even EFT corrections, and the more general mixed parity case.Both cases are captured by equation (120).

Generalization to Kerr
We can now generalize our conclusions to Kerr BHs of any spin.As before, we will prove that the corrections to QNM frequencies can be obtained from the parity even results, section IV B 3, by simply replacing λ (e) with |λ| = λ 2 (e) + λ 2 (o) .As in the previous section, the starting point will be equation (134).Following the procedure outlined in section IV B until equation ( 187), but keeping λ (o) ̸ = 0, we get: where we remind the reader that Φ 0 is the component of the corrected electromagnetic tensor we want to determine, ψ H is a Hertz potential for the background, λ = λ (e) +iλ (o) is a complex O(1) constant, P is the parity operator, and P the conjugate-parity-transform operator (see section II E).D and T 0 , defined in section IV B 3, are GHP operators with spin-weight s = 0.In the parity even section, we have proved that P commutes with all operators in (232): Equation ( 232) is almost in the same form as equation (210).However, in the place of A, B, differential operators that commute with complex conjugation, we have operators that commute with P.This motivates the definition of a new projection operator B δ , analogous to A δ defined in (212).For some complex function η, we have: As argued at the end of section II E, P preserves GHP weight, thus this is also true for B δ .Furthermore, note that all the properties listed for A δ at the beginning of section IV C, generalize to B δ , by replacing complex conjugation with P.
Because P commutes with all the GHP operators in equation ( 232), so does B δ .Acting with it on both equations, we get: where λ θ is the phase of λ.In both equations we see B δ Φ 0 , however, in the top equation it is related to B δ−λ θ ψ H , whereas in the bottom equation we relate it to B −δ ψ H .To get decoupled equations we must have B −δ ψ H ∼ B δ−λ θ ψ H .This condition is entirely analogous to what we explored in the toy model at the beginning of section IV C. We find two decoupled equations, analogous to (215): This reduces to the EOM for the parity even case, equation (192), under the transformations: 19 Thus, in agreement with the findings of section IV C 1, we prove that the corrections to Kerr QNM frequencies under mixed parity EFT corrections, can be obtained from the parity even case by taking λ (e) → λ 2 (e) + λ 2 (o) .We have proved that for any Kerr BH background, QNM frequencies are described by equation (120).

D. Results
In this section, we will solve the EOM numerically, to find the corrections to the QNM frequencies.In sec- 19 The factor of i 1∓1 2 is included because Φ ± 0 and ψ ± H are parity even / odd scalars, whereas B δ ± Φ 0 and B −δ ± ψ H are both parity even.This is purely technical and has no bearing on the final result.
tion (IV D 1) we will study the Schwarzschild limit of parity even EFT corrections, and then, in IV D 2 we generalize the result to Kerr BHs of arbitrary spin.As argued in IV C, at leading order in L/M , the results can be readily generalized to mixed parity corrections by replacing λ (e) with λ 2 (e) + λ 2 (o) 1. Schwarzschild limit FIG. 9.In this plot, setting λ (o) = 0, we show the corrected electromagnetic QNM frequency, for Schwarzschild BHs as a function of λ (e) .The plot lines show the solutions to equation (179) obtained non-perturbatively in L/M .The plot markers were extracted from figure 4 of [51].In that paper, the authors computed the QNMs using a third-order WKB approximation, valid in the limit ℓ ≫ 1. ℓ = 1 is outside this regime, thus, as expected, there is very poor agreement between our results for ℓ = 1.The situation improves for ℓ = 2 and ℓ = 3, but for large values of λ (e) the difference is still noticeable.As argued before, the EFT approximation is only valid in the limit L ≪ M , thus most of this plot is non-physical.
Using the direct pseudospectral methods outlined in the appendix A, we solved equation (179) for ℓ = 1, 2, m ∈ {−ℓ, • • • , ℓ}, by working non-perturbatively in L/M .We took ∼ 50 values of λ (e) ∈ [−1, 1] and obtained the corresponding QNM frequency for both polarizations, ω ± .As seen in section IV B 1, in the Eikonal limit, for Schwarzschild BHs, the imaginary part of ω 1 vanishes.Thus, for ℓ ̸ = ∞, we expect Im(ω ± 1 ) to be small.This implies that, for Schwarzschild BHs, Im(ω ± ) will be particularly close to Im(ω 0 ).Thus, in this regime, to appropriately estimate ω ± 1 , we must have very precise values of ω ± and ω 0 .In [51], the authors solved equation (179) for many values of λ (e) , using a third order WKB approximation.This approach is only valid in the large ℓ limit.Although several results in the literature argue that WKB methods produce reasonable results even for ℓ = O(1), the imprecision will necessarily grow the closer we get to ℓ = 1.In figure 9, we compare our results with theirs.Although we find rough agreement for ℓ = 2, 3, for ℓ = 1 we see a large difference.In fact, the difference between the two results at λ (e) = 0 is of the same order of magnitude as the maximum correction for λ (e) ∈ [−1, 1].This is not surprising given the discussion above.III.In this table, we compare three values for the Schwarzschild limit of ω1.The Direct column shows the correction estimated using equations ( 179) and (239), taking λ (e) = 10 −9 .The Perturbative column contains ω1, computed using the perturbative approach developed in section IV B 3, evaluated at a = 0.The last column has an estimation of this frequency deduced from the data in figure 4 of [51].The results in the first two columns were truncated to 10 significant digits, but we actually observed an absolute difference of at most ∼ 10 −19 .Furthermore, as expected, the corrections for the two polarizations are equal and opposite.As expected, due to the poor accuracy of the WKB approximation used in [51] at low ℓ, there is no agreement in the last column for ℓ = 1.The situation improves with increasing ℓ.

Parameters
We have proved, that at leading order in λ (e) , the EFT corrections to the QNM frequencies corresponding to each polarization should be equal and opposite.The results in figure 9 seem to be consistent with this notion.To accurately evaluate this, we must specialize to the small λ (e) regime.Assuming ω ± λ (e) is regular near λ (e) = 0, we have: (239) Plugging in λ (e) = 10 −9 we obtained very precise values for ω 0 and ω ± 1 .The background component describes an electromagnetic QNM of a Schwarzschild BH.This has recently been calculated with very high precision, using pseudospectral methods in [37].We find exact agreement between our results and the frequencies in table 1 of that paper.In table III, we show the results for ω ± 1 .We find that ω + 1 + ω − 1 ≲ 10 −19 .In the next section, working perturbatively in L/M , we will generalize the results to Kerr BHs.Setting a = 0, we get an independent result for ω ± 1 .Once again, we find agreement with at least 19 digits of precision.We can also estimate the value of ω ± 1 using the data in [51].We fit a second-degree polynomial to the data with small λ (e) , and take the derivative at λ (e) = 0.Because there isn't a lot of data in this region, we get a very rough estimate.This is included in the last column of III.As in the non-perturbative regime (see figure 9), the results are roughly consistent with ours for ℓ = 2, 3, but there is no agreement at ℓ = 1.

Kerr results
To solve equation ( 203), and obtain the QNM frequencies corresponding corrections, we employed pseudospectral methods, as described in the appendix A. The main results are summarized in figure 2 and table IV.We will discuss qualitatively these results below, but we will first establish their consistency and accuracy.In figure 4, we computed ω 1 for a BH with a = 0.999975M at varying grid sizes.We proved that, with increasing grid size, the numerical accuracy of ω 1 increased exponentially.This is in line with the expectation for pseudo-spectral methods.
One of our main consistency checks was done in the Schwarzschild limit.In section IV B 2, we derived the EOM for the QNMs in this regime, following an approach that is not perturbative in λ (e) .In table III, we compared ω 1 derived through that approach, with the solutions of equation ( 203), with a = 0. We found perfect agreement within the numerical accuracy of both approaches.Due to spherical symmetry, in this regime, QNMs and corresponding corrections are independent of m, with the degeneracy being broken for a > 0. This can be seen in figures 2 and 12.
On the other hand, our results are consistent with the Eikonal limit prediction of section IV B 1.
In figures 10 and 14 we show, accordingly, the ℓ = m and ℓ = −m QNM frequencies of Kerr BHs for a large range of ℓ values, together with the analytic large ℓ limit prediction.We find that, with increasing ℓ, the numerical curves converge to the analytic prediction.To quantify the convergence rate, in figures 11 and 15, we used a logarithmic derivative approach.We find that the real part of the background frequency converges at a rate of order ℓ −1 , whereas the imaginary part converges as ℓ −2 .This is consistent with the prediction in [14].Although we had no results regarding the expected convergence rate for ω 1 , we found that the real and imaginary parts approach the Eikonal limit as ℓ −1 .Qualitatively, our results are very similar to the ones obtained for the scalar case (see section III C).The crucial difference is the splitting of the QNM frequencies into two families, corresponding to the two possible polarizations of light.We have:

Parameters
We see that parity even and parity odd EFTs corrections change the QNM frequencies in the exact same way.Furthermore, we find that the corrections corresponding to the two polarizations are equal and opposite.This only happens because, in the RHS of equation (187), we were able to cancel all terms proportional to Φ i , leaving only terms involving Φi .Otherwise, we should expect the corrected frequencies to not be centered around ω 0 .This seems to be the case in [33][34][35].Presumably because the EFT corrects the background spacetime, the authors find that ω + 1 + ω − 1 ̸ = 0. Notwithstanding, the authors find that the mixing between parity even and parity odd QNMs is given by the square root of the sum, in agreement with (240).
Modes with ℓ = m, also dubbed slowly damped QNMs, have background frequencies, with the least negative imaginary part.This means they are the slowest decaying modes, thus dominating the spectrum of the BH ringdown.From these, for most of the parameter space, the ℓ = m = 1 is the lowest lying, However, at a ≈ 0.96M , Im(ω 0 ) for ℓ = m = 2 and ℓ = m = 3 cross the value for ℓ = m = 1, and these become the slowest decaying QNMs. 20.We did not see this tendency for modes with ℓ > 3, but it could be present sufficiently close to extremality.Coincidentally, except very close to extremality, for each value ℓ, ω 1 has the largest absolute real and imaginary part when ℓ = m.These hit a maximum for As expected, with increasing ℓ, the QNM frequencies approach the large ℓ limit line.Furthermore, in agreement with figures 2 and 12, in the near extremal limit, M ω0 → ℓ/2 and M ω1 → 0.
a ∼ 0.75M .Thus. these modes seem to be the most relevant if we wish to experimentally find an upper bound for λ.
In figure 12, we show the data of figure 2 in a parametric form, as a function of a/M .The plot markers represent QNM frequencies in the Schwarzschild limit, as obtained in section IV D 1, whereas the plot lines correspond to the QNM frequency when a ̸ = 0.There is a perfect agreement between the two.Particularly relevant, is the shape of the slowly damped QNMs.Their parametric lines cover a large region of the complex plane, varying particularly quickly in the near extremal lime.For these modes, should we exclude a ≳ 0.99M , we would lose approximately the last tenth of the parametric lines.Thus, in figure 13, we zoom in this region of parameter space, computing the QNMs for BHs with a/M ∈ (0, 99, 0.999975).As predicted in [38], as a → M , the background QNM frequencies seem to converge polynomially to ℓ/2.Analogously, in the same limit, we find that M ω 1 approaches 0 at a similar convergence rate. 21 Furthermore, we find that the absolute value of M Reω 1 grows with ℓ.This is consistent with the prediction of [38] for the background modes.Finally, note that at a ∼ 0.992, M Re ω 1 crosses 0.

V. DISCUSSION
In this paper, we computed the correction to BH QNMs, under scalar and electromagnetic, leading order EFT corrections.In the electromagnetic case, 21 In the ℓ = m = 1 and ℓ = m = 2 cases, modes seem to converge slightly quicker than expected, but, for ℓ = m = 3, the convergence rate is consistent with the expectation: ω ℓ=m=3 0 − ℓ/2 ∼ √ M − a.We think that, by getting even closer to extremality, we should get the correct convergence rate.It was harder to quantify the convergence rate of ω 1 as a logarithm derivative approach did not seem to converge for a ≲ 0.999975M .This is consistent with the results in the figure 10.We did not derive any theoretical results regarding the expected convergence speed of the EFT correction.Notwithstanding, we find that for both the real and imaginary parts of ω1, the convergence rate is consistent with ℓ −1 .Note that there is erratic behaviour of the EFT correction for low ℓ values, particularly for a/M = 3/7 and a/M = 4/7.This happens, because, as seen in figure 10, ω ℓ 1 crosses ω ∞ 1 , in this region, taking the logarithm of the difference to infinity.
we proved that the leading order parity even and parity odd EFT corrections have the same effect on QNM frequencies.Taking λ(e)/(o) to be the coupling parameter that control parity even / odd corrections, the combined effect of the two is controlled by λ 2 (e) + λ 2 (o) .Furthermore, we found that the degeneracy between the two polarizations of light is broken, splitting the QNMs into two families.At leading order, the corrections to the corresponding frequencies were shown to be equal and opposite.
Electromagnetic QNMs from a BH merger have not yet been detected.Thus, from an observational perspective, our results are not very relevant.However, as seen in [33][34][35], the main features generalize to pure gravity EFTs corrections.There, the degeneracy breaking between the two possible polarizations is particularly relevant.In theory, if we know the precise direction of a GW source, combining the signal from several GW detectors, we may disentangle the two GW tensor polarizations, see [68].
Thus, if both families of QNMs are excited in a BH merger, we may be able to disentangle the QNM spectrum into contributions corresponding to each polarization, and use the difference to infer an upper bound for the size of the EFT correction.
We used, scalar, and electromagnetic QNMs as toy models to learn how to, in general, compute the EFT corrections to gravitational QNMs of Kerr BHs.There, the challenge is that not only the EOM are corrected, but the background stationary metric is also corrected.[33][34][35] have overcome this challenge, by working perturbatively in the BH spin.In future work, for BHs with a/M ∈ (0, 0.999975).We find great agreement between the two results.Qualitatively, the frequencies share most features of scalar QNMs.In the Schwarzschild limit, modes with the same ℓ share the same frequency and correction.In the near-extremal limit, for ℓ = m modes, the background frequency approaches ℓ/2, whereas the correction approaches 0. we expect that it is possible to extend our method to compute the gravitational QNMs of BHs with arbitrary spin.The correction to the background metric can be determined numerically and then, the resulting EOM can be decoupled and solved with an extension of both our methods.Mapping the remaining factors to "diagonal tensors", the discretization of (A8) will take the form: where the Einstein summation convention was assumed.
When using pseudo-spectral methods, it is important to quantify the size of the numerical error.Because these methods converge exponentially with the grid size, we can simply solve the EOM using two different grid sizes.The absolute difference between the values obtained for each is a good estimate of the error size.In this paper, in most cases, we set the second grid to be roughly 10% bigger than the first.

Solving the background equation
Following the discretization procedure, we can now solve the EOM and find the QNMs.We will first describe how to solve the background equation, equation (A1a).As argued before, this equation separates into a system of two ODEs.Choosing a Chebyshev grid with size N x in the angular direction and a grid with size N r in the radial direction, the discrete representation of equation (A1a) is: A x (J, m, A ℓm , ω 0 ) • Θ = 0 , (A12a) A r (J, m, A ℓm , ω 0 ) • R = 0 . (A12b) Here A x and A r are matrix operators, J and m are parameters of the problem (see equations ( 14) and ( 46)), A ℓm and ω 0 are to be treated as eigenvalues, and R and Θ are the corresponding eigenvectors.Throughout the rest of this section, all underlined quantities are to be understood as vectors.The two equations depend on both eigenvalues, hence the system is not completely decoupled.Added difficulty comes from the fact that the second equation is quadratic in ω.
To alleviate this, we can start by solving the EOM in the Schwarzschild limit.As argued in section II C, in this limit, the angular equation can be solved analytically, yielding the spin-weighted spherical harmonics.A ℓm reduces to l(l + 1) − s(s + 1).The symmetry of the problem also removes the m dependence from the sec- In this image, we study how quickly do ℓ = −m QNM frequencies converge to the Eikonal limit prediction.As in figure 11, we find that the real and imaginary parts of the background converge as ℓ −1 , and ℓ −2 .This is consistent with the prediction in equation (151).Furthermore, we find that the real and imaginary of the correction converge to the theoretical prediction as ℓ −1 .
ond equation.Following section III.C of [19], we used a linearization approach to reduce (A12b) to a standard eigenvalue problem.The resulting discrete eigenvalue problem can be inserted in a numeric eigenvalue solver.We used Mathematica's Eigensystem method, selecting the Arnoldi method. 22To achieve accurate results, the use of Mathematica's arbitrary precision capabilities was crucial.In general, we used at least 200 digits of precision.
To tackle the case with arbitrary spin, we used the Newton-Raphson approach, (see section VI.A of [19]).This is an iterative method to find the roots of any map f : R d → R d .Define x 0 such that: Taylor expanding f (x) around x 0 , we get: where δx := (x − x 0 ) and δf (x) := ∂f (x) ∂x .By definition f (x 0 ) = 0, thus, we can approximate δx with: We update x → x − δx, and find a new estimate of δx, iterating the process until convergence.As long as the initial seed is close enough to x 0 , this method has quadratic convergence.However, if the seed is too distant from x 0 , the method will either not converge, or converge to a different root of f .
To apply this method to (A12), we must first fix some leftover gauge freedom.The vectors Θ and R can be rescaled by constant factors, thus, we must supplement the EOM with two normalization conditions.We choose two vectors n (x) and n (r) , and require: Taking x = (Θ, R, ω 0 , A ℓm ), we use equations (A12) and (A16) to define f .The Newton-Raphson update equation will be: Now, to find the background QNM frequency for a background of any spin, we just need an accurate enough initial guess.We take the Schwarzschild limit frequency and use it as a seed for a very slowly rotating BH.After convergence, we can use the resulting frequency to obtain the QNM of a BH with a slightly larger spin.We iterate this, increasing J gradually, until we cover the whole range of Kerr BHs, J ∈ [0, 1).
Deciding how much to increase J at each iterative step, is subtle.To aid in this decision, and speed-up convergence, we can estimate dx 0 /dJ.Including J in equation (A13), we have f (x 0 , J) = 0. Thus, we can treat x 0 as a function of J, parametrizing the level curves of f .Along these curves, we have: Solving with respect to dx 0 /dJ we get: Plugging this in the definition of the derivative, we can estimate x 0 (J + δJ): x 0 (J + δJ) ≃ −x 0 (J) − (δf ) −1 • ∂f ∂J δJ . (A20) In our calculations, we fixed the size of δIm(ω 0 ), deducing the corresponding size of δJ.

Solving the perturbation equation
Using the background QNMs, we can now solve the perturbation equation, equation (A1b).Following the prescription in section A 1 we obtain the corresponding discrete version.Because the EOM is a 2 dimensional PDE, Φ (0/1) are represented by N x × N r matrices, and L (0) , L (1) and δL will be represented by 4-tensors with size N x × N r × N x × N r .Fixing normalization freedom as was done (A16), we get: where we assumed the Einstein summation convention.These tensors depend on the same parameters as the background and on (ω 0 , A ℓm ).Furthermore, the separability condition at 0 order translates to Φ (0) = Θ⊗R.Using equation (A16), the normalization condition (A21b) simplifies to n (x) • Φ (1) • n (r) = 0. Using Mathematica's Flattening operator, we get a vector representation of Φ (0/1) .This in turn implies a matrix representation of the 4-tensors above.In this basis, the discrete EOM take the form of a linear system: where: x = Φ (1) , ω 1 , δL (1) • Φ (0) T n (r) ⊗ n

(A23)
Underlined quantities are N x • N r vectors whereas not underlined quantities are scalars or (N x • N r ) × (N x • N r ) matrices, depending on the context.The system (A22) can then directly be solved with Mathematica's LinearSolve method.In general M is non-singular, thus, we expect that there will be a single ω 1 corresponding to each background frequency ω 0 .
Chebyshev differentiation matrices are very dense, thus, in general, the discrete systems we obtain have very high condition number.Thus, to solve the EOM we must work with very high-precision arithmetic.This is particularly relevant in the near extremal limit, where the condition number increases very fast.M and b are functions of ω 0 and A ℓm , thus the system's precision is determined by the precision of these eigenvalues.To accurately solve the EOM in the near extremal regime, we needed these eigenvalues to have at least 70 digits of precision.To achieve such feat, we had to solve the background EOM using very large Chebyshev grids.In the most extreme scenario, we used N x = 250 and N r = 650.When solving the background EOM, this implies that each step of the Newton Raphson method involves solving a ∼ 900 × 900 system of equations, which is feasible with current computing capabilities.However, when discretizing the first order PDE, rank of M will be N r • N x ∼ 160000.Even though this matrix is not very dense, solving this system would be very challenging.To work around this issue, we can use a clever form of interpolation to reduce the rank of M. As mentioned in section A 1, there is a correspondence a Chebyshev collocation grid and the spectral decomposition in a Chebyshev polynomial basis.In fact, these coefficients are related by a linear map that can be found in [21].To project the perturbation equations into a smaller Chebyshev grid, we find their spectral representation, truncate the spectral sum at a lower order, and then, map it back to a Chebyshev collocation grid.

FIG. 1 .
FIG.1.This image summarizes the EFT corrections to low ℓ 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 ≲ 0.999M in most cases, or a ≲ 0.999975M for ℓ = m ̸ = 0 modes.We found that, in general, the imaginary part of the background frequency increases with ℓ.Thus, the lowest lying QNMs should have ℓ → ∞.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 ℓ = m ̸ = 0 modes M ω0 → ℓ/2, as predicted in[38].This feature extends to the EFT correction: ω1 → 0 as a → M .

FIG. 4 .
FIG.4.In this plot we show that our numerical approach converges exponentially with increasing grid size.Here, Nx, Nr 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 ℓ = m = 1, for multiple values of Nx and Nr.Then, we compared them with a reference ω1 obtained with the largest grid, plotting the difference.On the left panel, we fixed Nr = 250 and varied Nx, whereas on the right panel we fixed Nx = 75, and varied Nr.The electromagnetic case was particularly difficult to compute, to achieve accurate results we had to obtain ω0, with 70 digits of precision.

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

FIG. 11 .
FIG. 11.In this image we study how quickly the QNM frequencies approach the Eikonal limit.On the y-axis, we see the logarithmic derivative, γ = d log ω ℓ − ω ∞ /d log(ℓ).Panels a) and b), show the convergence of the real and imaginary part of the background frequency, whereas c) and d) have the real and imaginary part of the EFT correction.Because we expect Reω to grow linearly with ℓ, we normalized the difference diving by ℓ.We have M ω ℓ − ω ∞ ∼ ℓ −γ .From equation (151) we expect that the real part of the background frequency has γ → −1, while the imaginary part should have γ → −2.This is consistent with the results in the figure10.We did not derive any theoretical results regarding the expected convergence speed of the EFT correction.Notwithstanding, we find that for both the real and imaginary parts of ω1, the convergence rate is consistent with ℓ −1 .Note that there is erratic behaviour of the EFT correction for low ℓ values, particularly for a/M = 3/7 and a/M = 4/7.This happens, because, as seen in figure10, ω ℓ1 crosses ω ∞ 1 , in this region, taking the logarithm of the difference to infinity.

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

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

FIG. 14 .
FIG. 14.Here we show the QNM frequencies of ℓ = −m electromagnetic modes for several values of ℓ.In panels a) and b), we have the real and imaginary part of the background frequency, with panels c) and d) showing the corresponding correction.The thick dark line shows the Eikonal limit prediction.As expected, the QNM frequencies converge to the asymptotic prediction.
FIG.15.In this image, we study how quickly do ℓ = −m QNM frequencies converge to the Eikonal limit prediction.As in figure11, we find that the real and imaginary parts of the background converge as ℓ −1 , and ℓ −2 .This is consistent with the prediction in equation (151).Furthermore, we find that the real and imaginary of the correction converge to the theoretical prediction as ℓ −1 .

TABLE IV .
Numerical values for the QNM frequencies and corresponding EFT correction of electromagnetic modes.All the included digits are significant.