Analytically Separating the Source of the Teukolsky Equation

Recent gravitational wave detections from black hole mergers have underscored the critical role black hole perturbation theory and the Teukolsky equation play in understanding the behaviour of black holes. The separable nature of the Teukolsky equation has long been leveraged to study the vacuum linear Teukolsky equation; however, as theory and measurements advance, solving the sourced Teukolsky equation is becoming a frontier of research. In particular, second-order calculations, such as in quasi-normal mode and self-force problems, have extended sources. This paper presents a novel method for analytically separating the Teukolsky equation's source, aimed to improve efficiency. Separating the source is a non-trivial problem due to the angular and radial mixing of generic quantities in Kerr spacetime. We provide a proof-of-concept demonstration of our method and show that it is accurate, separating the Teukolsky source produced by the stress-energy tensor of an ideal gas cloud surrounding a Kerr black hole. The detailed application of our method is provided in an accompanying \textit{Mathematica} notebook. Our approach opens up a new avenue for accurate black hole perturbation theory calculations with sources in both the time and frequency domain.


INTRODUCTION
Since the first detection of gravitational waves from a black hole binary [1] in 2015, gravitational wave astronomy has been a rapidly growing field.The current catalogue of binary detections [2] is centred on binaries where the two compact objects have a similar mass.However, recently, the LIGO/Virgo/KAGRA collaboration detected a binary with mass-ratio of ε = µ M ∼ 0.04 [3], where M and µ are the masses of the primary the secondary binary component respectively.Gravitational astronomy is now not confined to the comparable mass regime; the so-called intermediate-mass-ratio inspirals (IMRIs, with mass ratio 10 −2 ≳ ε ≳ 10 −4 ) are increasingly more detectable [4,5].Additionally, millihertz frequency gravitational waves will be detectable with space-based gravitational wave interferometers, sensitive to lower mass ratios than groundbased detectors.In the 2030s, three space-based gravitational wave interferometer missions are planned: LISA [6], TianQin [7], and Taiji [8].Extreme-mass-ratio inspirals (EMRIs, 10 −5 ≳ ε ≳ 10 −8 ) are a key target for these missions [9].
Numerical Relativity (NR) [10] has had great success in accurately modelling comparable mass binaries.Progress is being made to extend the catalogue of NR waveforms into the IMRI regime [11,12].Synergies with black hole perturbation theory can also improve efficiency [13,14].However, these are challenging tasks due to the disparate length scales of such systems.In the EMRI regime, it is unrealistic to expect such methods to be capable of filling the parameter space of initial conditions in the LISA mission time frame.
Alternatively, black hole perturbation theory provides an ideal method for approximating the binary spacetime in the small mass-ratio limit.The evolution of the binary can be modelled using the self-force approach [15].In Refs.[16][17][18], a two-timescale approximation [19,20] was implemented for a self-force evolution of a binary, including (the dissipative) second-perturbative-order self-force.This model was limited to a binary of two Schwarzschild black holes in a quasi-circular inspiral.The waveforms were produced for a range of mass ratios and showed outstanding agreement with NR (before the inner-most stable circular orbit), even in the ε ∼ 1  10 regime.These results show that self-force waveforms have incredible potential for gravitational wave science.Crucially, the dissipative piece of the second-order self-force is a necessary contribution for high accuracy waveforms.There is also growing interest in second-order quasi-normal mode calculations [21][22][23][24][25][26][27][28][29].
The recent second-order self-force results are limited to a Schwarzschild background black hole, whereas the primary black hole of astrophysical EMRIs is expected to have a significant spin (a ∼ 0.9M ).Adding the linear in-spin contributions perturbatively [30,31] is straightforward.However, including the full nonlinear spin of the primary object is one of the outstanding problems in second-order self-force.Working non-linearly in spin involves the background spacetime being a Kerr black hole [32].The Kerr metric [32], in Boyer-Lindquist coordinates [33], is where M is the black hole mass, a is the angular momentum per unit mass, Σ = r2 + a 2 cos 2 [θ], and ∆ = r 2 − 2M r + a 2 .Kerr spacetime is stationary and axially symmetric, but not spherically symmetric; hence, there is radial (r) and polar angle (θ) mixing in the metric from (1).This lack of symmetry makes the linearised Einstein field equations generally non-separable in Kerr; attempting to reduce the linearised Einstein field equations, a coupled set of PDEs, to a coupled set of ODEs by separation of variables seems fruitless.However, a separable field equation does exist, the Teukolsky equation [34,35].Solving the sourced Teukolsky equation separably requires separating the angular and radial dependency in the source.
In this paper, we produce a formalism for analytically separating general sources of the Teukolsky equation.Our method will be applied to produce accurate second-order self-force and quasi-normal mode calculations in follow-up papers.We expect our method will provide key efficiency savings in the second-order self-force problem, where a substantial amount of mode coupling is present.Before presenting our formalism, we briefly introduce black hole perturbation theory and the separability of the Teukolsky equation in the remainder of this section.In section II, we present our method for decomposing the source (for a summary, see Fig. 3).In section III, we apply our method to a toy example to show it accurately decomposes sources.Finally, in section IV, we make our concluding remarks.

A. Black hole perturbation theory and separability
Black hole perturbation theory involves solving the Einstein field equations perturbatively around a background black hole spacetime.We write Einstein field equations in natural units as where G ab is the Einstein tensor, g ab is the metric1 , and T ab is the stress-energy tensor.The metric and stressenergy tensor are expanded in orders of the small parameter ε, We take g ab , the background metric, to be the Kerr metric.Therefore, T (0) ab = 0, as Kerr is a vacuum solution.Under the expansions in Eqs. ( 3) and (4), the Einstein field equation, Eq. ( 2), can be expressed as a hierarchical set of linear field equations in ascending order of ε [17,36], δG ab [h ab ] = 8πT (2) ab , h ab ], δG ab [h ab ] = . . .
where δG ab is the linearised Einstein operator and δ 2 G ab is the quadratic Einstein operator [36].
It is not immediately obvious that Eq. ( 8) is separable and relates to Kerr spacetime being Petrov type-D [40] and vacuum 2 .To show Eq. ( 8) is separable, one can choose a tetrad and convert to master Teukolsky equation form [34,35,39].In the Kinnersley tetrad [42], the master Teukolsky equations are ab ], where Ôs is the spin s master Teukolsky operator and The separability ansatz relies on the time dependence of the Teukolsky master variable, ψ (1) 0 4 , being expressed as proportional to e −iωt (where ω may be complex); for example, working in the frequency domain.The spin s master Teukolsky operator can then be expressed as where [39], in Boyer-Lindquist coordinates, where K := (r 2 + a 2 )ω − am, and where χ := cos[θ] and A := s λ ℓm + 2amω − a 2 ω 2 (and the eigenvalues s λ ℓm depend on aω).The eigenfunctions of Θ s are spin-weighted spheroidal harmonics3 [34,35], More generally, the Teukolsky equation can be shown to be separable in any principle-null direction aligned tetrad in Kerr without using coordinates [39,43].The Teukolsky equation can be written as where R and S are operators that commute; hence, the equation is separable.
Research has largely focused on solving the first-order vacuum Teukolsky equation [44][45][46] or point-particle source problem [15,39].Hence, previous research has generally not required efficient separation of the right-hand side of the Teukolsky equation into modes of the separable Teukolsky equation.As calculations move to second order, sources are unavoidable.The second-order Teukolsky equation comes in two forms, the Campanelli-Lousto second-order Teukolsky equation [47] and the reducedsecond-order Teukolsky equation [36,48].Both equations are sourced by the stress-energy perturbations and a quadratic operator acting on the first-order metric perturbation.For example, the reduced second-order Teukolsky equation takes the form 4L ] = S 4 8πT (2) ab , h ab ] , where ψ (2) ab ].Note that ψ 0L and ψ (2) 4L are the linear in h (2) ab part of the full second-order Weyl scalars, where ψ ab , h ab ] and ψ ab , h ab ] are the quadratic in h respectively [36,47]; the operators δ 2 ψ 0 and δ 2 ψ 4 are given in Ref. [49].The ab ] parts of Eqs. ( 17) and ( 18) make the source non-compact, extending from the horizon of the primary to future null infinity.
We need a precise definition for separating the source of the Teukolsky equation.Take the master Teukolsky equation with a generic source f [t, r, θ, ϕ], for the spin-s master Teukolsky variable ψ.In order to solve the equation separably, ψ and f must be expressed in a mode decomposition form: Then, Eq. ( 21) reduces to for each ℓ and m.
As spin-weighted spheroidal-harmonics are an orthonormal set [50] of eigenvectors on S 2 , one can separate the source of the Teukolsky equation of spin weight s, using where s S * ℓm is the bi-orthogonal dual [51] of the spinweighted spheroidal harmonic which simplifies to the complex conjugate when ω is real.f ℓm [t, r] in Eq. ( 25) will then satisfy Eq. ( 23).As spin-weighted spheroidal harmonics generally have no known closed form, integrating Eq. ( 25) is challenging.One could integrate numerically, but f [t, r, θ, ϕ] is not generally separable in the radial and polar angle coordinates.Therefore, the integral will need to be computed numerically at each radial point on a grid for each ℓ and m mode.This method has been used to separate the source of the Teukolsky equation in Ref. [29] for a quadratic quasi-normal mode calculation.However, the inefficiency of the numerical integrals may be problematic for self-force calculations where multiple modes must be calculated, and second-order calculations must fill a four-dimensional parameter space of initial data to create the waveform templates for LISA data analysis.
Alternatively, one could avoid separating the source of the Teukolsky equation by solving the Teukolsky equation in PDE form [52,53].However, numerical methods for solving PDEs convert the PDE into a system of ODEs.Hence, it is likely that leveraging the separability of the Teukolsky equation will develop the most efficient algorithm for separating the Teukolsky equation into ODEs.

II. DECOMPOSING GENERAL FUNCTIONS IN KERR SPACETIME
Our formalism is designed to separate general sources in Kerr into the sum of the product of radial functions and spin-weighted spheroidal harmonics.Essentially, we separate the radial and polar angle dependency 4 and express the angular dependence in terms of spin-weighted spheroidal harmonics.
An immediate simplification we can make is reexpanding the spheroidal harmonics into spin-weighted spherical harmonics.Spin-weighted spherical harmonics are an easier basis of functions to work with because they are closed-form and have simple spin-weight raising and lowering operators associated with them [54] 5 .Leavers method [56] is the most accurate method of expressing spin-weighted spheroidal harmonics but does not expand them in terms of spin-weighted spherical harmonics directly.Press and Teukolsky [57] were the first to expand spin-weighted spheroidal harmonics into spin-weighted spherical harmonics.Ref. [58] then found a more efficient method, which has been implemented in the Black Hole Perturbation Toolkit [59].The expansion in terms of spherical harmonics converges to Leaver's method solutions and is invertible.But, to work with either expansion, one must truncate at a finite order.Hence, for the remainder of this paper, we assume any inputs that are generally expressed in terms of spin-weighted spheroidal harmonics can be re-expanded in terms of spin-weighted spherical harmonics.Similarly, it is sufficient to express the source in terms of spin-weighted spherical harmonics as we can re-expand it in terms of spin-weighted spheroidal harmonics.That is, we can simplify our goal to expressing the source as where s Y ℓm [θ, ϕ] are the spin-weighted spherical harmonics.The rest of this section is dedicated to expressing the source of the Teukolsky equation in the form of Eq. ( 26).To simplify this task, we use the Newman-Penrose [60], Geroch-Held-Penrose (GHP) [61], and Held [62] formalisms.We also find there is a preferred tetrad and coordinate scheme to express the source in Eq. ( 26) form.In such a tetrad and coordinates, almost all the Kerr background quantities are naturally expressed as single spin-weighted spherical harmonics, and all angular derivatives become spin-weight raising/lowering operators.One problematic background quantity is unavoidable, ρ (and ρ), where As the radial and angular dependence appears on the denominator, ρ does not naturally separate into a radial and angular function.In Sec.II E, we separate the radial and angular dependence in Eq. ( 27) using a Fourier expansion.Combining these methods, we produce a complete formalism for expressing the source as a convergent sum of spin-weighted spherical harmonics (Eq.( 26)).
A. Spin-weighted spherical harmonics The spin-weighted spherical harmonics were defined in Ref. [54], and we follow their conventions up to a minus sign in the definition of the spin raising and lowering operators (similarly to Ref. [38]), when acting on a spin-s weighted quantity.One can define the spin-s spherical harmonics from the scalar harmonics, The following relations then hold [54], The product of two spin-weighted spherical harmonics can be expressed as the sum of spin-weighted spherical harmonics, where C ℓ,m,s ℓ1,m1,s1,ℓ2,m2,s2 is equivalent to a surface integral, As spin-weighted spherical harmonics are related to Wigner-D matrices, Eq. ( 36) can be evaluated algebraically using 3j symbols [38,63,64],

B. The Newman-Penrose formalism
The Newman-Penrose (NP) formalism [60] utilises a basis of four null vectors, called a tetrad, to express curvature quantities.The NP basis vectors are labelled where indices in square brackets are tetrad indices.The vectors l a and n a are real, and m a is complex.Overbars denote a complex conjugation.Conventionally, for a positive metric signature, the orthonormal relationship of the tetrad takes the form As the tetrad is orthonormal, following Eq.( 39), the metric can be expressed as [60], The NP formalism uses Ricci rotation coefficients to express the connection [65], There are 24 independent components of Ricci rotation coefficients.In the NP formalism, the components are expressed as 12 complex scalars, known as spin coefficients: In vacuum spacetimes, the Ricci curvature is zero.The Weyl tensor contains the vacuum curvature.In the NP formalism, the ten degrees of freedom of the Weyl tensor are expressed using five complex scalars, known as Weyl scalars: Kerr spacetime has two principal null vectors (two pairs of principal null vectors which coincide) [40,66].That is, Kerr is Petrov type D. A tetrad can be chosen such that l a and n a are tangent to the principal null directions and four of the Weyl scalars and four spin coefficients then vanish [66,67], These simplifications were used to help derive the Teukolsky equation [35].
The covariant derivative is also expressed using tetrad components in the NP formalism:

C. The GHP formalism
The GHP formalism builds on the NP formalism, helping to represent the symmetry in principal null direction aligned tetrads in Petrov type D spacetimes (such as the Kinnersley [42], Carter [68], and Hartle-Hawking [69] tetrads).While constraining l a and n a to point in the same direction, two tetrad rotation degrees of freedom remain unconstrained.These freedoms can be associated with spin and boost transformations and are isomorphic to the group of multiplication by a complex number, ϑ [70].To express the weight of a GHP quantity f we use the notation f ⊜ {p, q}, meaning that under a spin and boost transformation p and q are known as GHP weights that can be equated to the spin (s = 1 2 (p − q)) and boost (b = 1 2 (p + q)) weights.Products of two tensors with weights a, b and c, d produces a tensor of weight a + c, b + d.On the other hand, the addition of two quantities can only be performed if they are of the same weight.
In Petrov type-D principle null direction aligned tetrads, there is a freedom to interchange l a and n a .GHP introduced a prime operation to represent the interchange l a → n a , n a → l a , m a → ma , and ma → m a .The prime operation affects the GHP weights accordingly, f ′ ⊜ {−p, −q}.Complex conjugation also affects the GHP weights: f ⊜ {q, p}.
In the GHP formalism, half of the NP spin coefficients are relabelled using the prime operation notation 6 , The GHP weights of the spin coefficients (and their primes) follow directly from the weights of the tetrad vectors (using Eqs.(42)) 7 ; that is, The spin coefficients ϵ, ϵ ′ , β, and β ′ do not have welldefined weights.Similarly, the NP derivative operators do not have well-defined weights.GHP found by combining these poorly defined weight quantities, one can produce derivative operators with well-defined weights, where η ⊜ {p, q}.Equation ( 54) respectively have boost and spin weights of Two clear advantages of the GHP formalism are that the equations are more condensed than in NP form, and they offer a straightforward consistency check by checking that the weights of an equation are consistent.
We will also make use of the GHP commutation relations [61] (with Eqs. ( 48) and ( 49) imposed), We apply the commutation relations, Eqs. ( 55) to (57) (and their complex conjugates) to commute all ð and ð ′ derivatives to the right, and all Þ ′ derivatives to the right of Þ.This simplifies equations, and this ordering will place all spin-raising and lowering operators to the right of radial derivatives in section II D 2.

D. The Held Formalism
The Held formalism builds on the GHP formalism, producing an algorithmic integration technique for a section of the GHP equations (the Ricci and Bianchi identities in GHP form).The Held integration method has recently been applied to solve black hole perturbation theory problems [48,70,71].Here, we review the main aspects of the Held formalism and its application in the Kinnersly tetrad [42].In Boyer-Lindquist coordinates, the Kinnersley tetrad is I also comment on the advantages of the Held formalism and when the integration technique can be used.
The Held formalism leverages the Ricci identities (called "Field equations" in [62]), Bianchi identities, and commutation relations to build a set of equations for derivatives of background Kerr quantities.For example, in a vacuum Petrov type-D spacetime, the R [4] Ricci identity [66] gives (62) Equation ( 62) can be interpreted as the integral of Þρ.
Using this identity to solve differential equations involving ρ and Þ is known as Held integration [62].For example, one may solve the complex ordinary differential equation.

ÞA[ρ]
where A[ρ] is an unknown function of ρ and B • is a function independent of ρ (ÞB • = 0, all variables labelled with • superscripts share this property and will be referred to as Held scalars8 ).Using Eq. ( 62) one can find the solution, where C • is any function independent of ρ.However, Held integration does not apply to general GHP equations.When an alternative derivative operator (such as Þ ′ , ð, or ð ′ ) acts on an unknown variable in a differential equation, the Held method of integration is not helpful.For example, cannot be solved using Held integration.Within these limitations, the Held integration method is useful for solving certain problems in black hole perturbation theory [48,70,71].
To extract the ρ dependency from the other spin coefficients, Held integrated various Ricci identities, Bianchi identities, and commutation relations [62], finding which defines the Held scalars There is an additional Held scalar that can be derived from the spin coefficient ρ, with At this point, we will deviate from the conventional Held formalism notation.In the Kinnersley tetrad we can express the Held quantities in a coordinated form.In Boyer-Lindquist coordinates ({t, r, θ, ϕ}) [48,70], Clearly, ρ ′• and Ψ • are coordinate invariant.In the following, we replace ρ ′• and Ψ • with −1/2 and M respectively.Note, this relabelling does not have a well-defined spin and boost weight, but this does not affect the form of 9 Held integration is still limited to equations where only Þ acts on the unknown quantity, whether using Held operators or GHP operators; like Eq. ( 65), ðD[ρ] = E • ρ 3 cannot be solved for D[ρ] using Held integration.
Bianchi identities [62]).Held used the remaining Ricci and Bianchi identities (and commutation relations) to derive simplifications for the Held derivatives acting on the Kerr-Held scalars, Þ′ 1. Held derivatives and spin-raising/lowering operators A further convenience of the Held formalism is that ð and ð′ relate to the spin-raising and lowering operators of the spin-weighted spherical harmonics.Next, we comment on how this relationship is not unique to the Held formalism and how the Kinnersley tetrad is one of a select class of tetrads where the GHP derivatives ð and ð ′ relate to the spin-raising and lowering operators.
In the Kinnersley tetrad, it is not challenging to extract the spin-raising and lowering operators from the GHP derivatives ð and ð ′ (Eq.( 54)) when expressed in coordinate form.First, express the GHP derivatives ð and ð ′ in NP form using Eq. ( 54).From Eq. ( 54), we see the relevant NP quantities are δ, β, α, and their complex conjugates.In the Kinnersley tetrad and Boyer-Lindquist coordinates, Eqs. ( 58) to (61) give Noting the common factor of 1 √ 2 ζ in Eqs. ( 81) and (82), and using Eq. ( 54), it is straightforward to show the GHP ð relates to the spin raising operator ð, The Kinnersley tetrad is not the only tetrad where GHP ð and ð ′ relate to spin-raising and lowering operators.This can be seen by examining the class III tetrad rotations [66] (under which the direction of the vectors l a and n a are unchanged and, therefore, remain aligned with the principle null directions), where A and θ are real functions.The relevant spin coefficients (see Eq. ( 54)) transform as [66], and the relevant spin coefficients in the Carter tetrad are Examining Eqs. ( 87) and ( 88), in the Carter tetrad, the cot[θ] terms in β and β ′ do not have the same coefficient as the θ derivative in δ.Hence the cot[θ] terms are not eliminated by the conversion to a spin-raising operator.Also, the a sin[θ] terms in Eq. ( 88) will remain after expressing the θ derivatives in terms of raising and lowering operators.cot[θ] and 1 sin[θ] are singular at the poles, and when expressed as a sum of spin-weighted spherical harmonics, the sum does not converge.Therefore, there are preferred tetrads for expressing Kerr quantities such that they are smooth on the 2-sphere: the Kinnersley tetrad and tetrads related to the Kinnersley tetrad by Eq. ( 84) with δθ = δA = 0.
While the Kinnersley tetrad in Boyer-Lindquist coordinates avoids singularities in ð and ð ′ , we can see from Eqs. ( 54) and ( 58) that Þ contains 1 ∆ terms.These terms are singular at the horizons as ∆ = 0 for the outer and inner horizons, respectively.To avoid singularities at the horizon and make Þ and Þ′ as simple as possible, we work in Kerr-Newman coordinates [48].The outgoing Kerr-Newman coordinates, u, r, θ, ϕ * , are related to Boyer-Lindquist coordinates as follows, [48,74] In these coordinates [48,70], Using Held's definitions [62] in the Kinnersley tetrad, it has been shown in Ref. [72], and The null vectors m a and ma in Eqs. ( 94) and ( 95) take the same form as in Boyer-Lindquist coordinates (as do the spin-coefficients as they are scalars which depend on only r and θ).The Held derivatives ð and ð′ then relate similarly to the spin raising and lowering operators, with the subtlety that the spin raising and lowering operators now act on The ∂ u derivatives are trivial as the background is timeindependent and we assume the perturbations have a time dependency proportional to e −iωt (converting the time coordinate dependence is straightforward, e −iωt = e −iω(u+r * ) ). Eqs. ( 98) and ( 99) are similar to Chandrasekhar's operators −s L ω † , +s L ω up to a trivial normalization [66].

Converting to coordinates and modedecomposition
To convert expressions in the Held formalism into Kerr-Newmann coordinates ({u, r, θ, ϕ * }), we first express the derivatives Þ and Þ′ using Eqs.( 96) and (97).The resulting expression presents all the ρ and ρ dependence explicitly 10 .Next, we deal with the remaining angular dependence and ð and ð′ operators.
We assume the inputs (A) are expressed as a convergent sum of spin-weighted spherical harmonics with coefficients that can depend on r, ρ, and ρ: where A ℓm [r, ρ, ρ] are polynomials with positive and finite powers of ρ and ρ; that is, for some finite i and j.Note that ϕ * derivatives, like time derivatives, are trivial as s Y ℓm ∼ e imϕ * ; that is, While we can use Eq.Eqs. ( 75) and ( 76) to express the Held scalars in coordinate form, for a mode decomposition, it is more useful to express them in terms of spin-weighted spherical harmonics 11 , where we have used with The resulting source expression will have products of many spin-weighted spherical harmonics.To combine the spin-weighted spherical harmonics, we use Eq.(35).The resulting expression for the source will be of the form Where fℓm [r, ρ, ρ] is a polynomial with positive and finite powers of ρ and ρ; that is, for some finite i and j.Next, we must express ρ and ρ in terms of spin-weighted spherical harmonics.

E. ρ expansion
So far, our method does not completely express quantities separably as factors where ρ and ρ are present; see Eq. (27).The mixing between the radial and polar angle coordinate in the denominator of Eq. ( 27) makes separating the r and θ dependency non-trivial.In this subsection, we show the r and θ dependency in ρ i ρk (for any positive integer i and k) can be separated using a Fourier expansion 12 .
As Eq. ( 27) is a trigonometric function, a Fourier expansion seems like a natural approach.Eq. ( 27) is periodic in θ (ρ[r, θ] = ρ[r, θ + 2nπ] for any r and θ, where n is an integer).Also, Eq. ( 27) is complex and non-singular for r ≥ r − (where r − is the radius of the inner horizon).To assess the real and imaginary behaviour of Eq. ( 27), one can write as All factors of ρ i ρk for i ̸ = k can also be written in a manner where the denominator is real: and for k ≤ i, This form is useful as the Fourier expansion of the numerators in Eqs. ( 106), ( 108) and (109) are trivial.Hence, only the Fourier expansion of the denominator, which is real, is required.A common denominator factor of Σ is desirable to minimise the number of Fourier expansions required.Using Eq. ( 106), we can express Eq.(105) as for some i, which has collected all of the factors of ρ and ρ outside of the sum in the form of 1 Σ i .Therefore, we require only one Fourier expansion, that of 1 Σ i .For example, take the function To extract the same ρρ dependency, we note that the highest power of either ρ or ρ is ρ4 .Hence, we write this function as, Using Eq. ( 27) to replace the inverse powers of ρ and ρ gives We can use cos[θ] = 4π 3 0 Y 10 and Eqs. ( 35) and ( 37) to combine the angular dependence with that in A • , B • , and C • , allowing us to express f in the form of Eq.(110).

Fourier series
A Fourier series expresses periodic functions using the orthogonality of the trigonometric functions sin[nθ] and cos [nθ].Truncated to order k, over an interval of 2π, the Fourier series of a function A[θ] can be expressed as where The quantity 1 Σ i is an even periodic function over [−π, π]; hence, b n = 0. Also, as Σ i is a function of only even powers of cos[θ], a 2n+1 = 0 for any positive integer n.Note that 1 Σ i is not singular in the domain r ∈ [r + , ∞) so we expect our Fourier series to converge.
We now assess the accuracy of a truncated Fourier series.For example, take the expansion of Σ −6 to order k = 8,  114)) is over-plotted as a dashed yellow line.Bottom: the fractional disagreement between the two plots is shown.There is less than 0.5% error, and the error decreases for increasing radius (see Fig. 2).calculated using Mathematica [78].Equation (118) accurately approximates Σ −6 , as can be seen in Fig. 1.For a = 0.9, r = 1.43M , and k = 8, the error is < 0.5%.For larger r, the approximation becomes more accurate (similarly to the blue line in Fig. 2).The most inaccurate region is for small r and near the poles.For a = 0.9, the outer black hole horizon is r + ≈ 1.436M (Eq.( 89)).For higher values of a, the truncated Fourier series is less accurate near the horizon 13 .Considering the accuracy requirements of LISA and current estimates on the spins of supermassive black holes, the expansion used in Fig. 1  F [Σ −6 ,k=12] − 1|, for different digits of precision used (black line is 8, red line is 32, and blue line is 128 digits of precision), for a = 0.9 and θ = π 4 .The plots show a clear dependence of the error on the number of digits of precision used for large r.This problem does not occur for k < 2i when approximating Σ −i . is probably sufficient.
To achieve arbitrarily high accuracy, one can increase the number of terms in the Fourier expansion, k.However, when attempting to approximate Σ n with a Fourier series with k ≥ 2n, errors can occur for large r if sufficiently high precision is not used, as shown in Fig. 2. We examine the simplest example, Σ −1 , to understand why large errors are encountered when using low precision.The Fourier series of Σ −1 with k = 2 is Naively, the Fourier series appears to scale as ∼ r 0 for large r, whereas 1 Σ ∼ r −2 .However, r − √ a 2 + r 2 → 0 for large r, so this contribution is suppressed, so the Fourier series scales as r −2 , as expected.However, resolving the suppression of such terms numerically requires high precision.With insufficient precision, an error ∝ r 2 is incurred.
In practice, using such high precision will slow calculations down.Nonetheless, for EMRI models for LISA data analysis, the accuracy requirements are likely sufficiently modest that Fourier expansions with k < 2i can be used; that is, high precision evaluation will likely not be required.

F. Combining the spin-weighted spherical harmonics
Finally, we can replace ρ i ρi = Σ −i factor in Eq. ( 110) with the Fourier expansion truncated to some finite order.
To convert the trigonometric quantities in the Fourier series with spin-weighted spherical harmonics, we can use Eqs.( 103) and (104).By again using Eqs.(35) and (37), we can combine the spin-weighted spherical harmonics to express the source as To solve the separable master Teukolsky equation, Eq. ( 24), we need to express the source as a spin-weighted spheroidal expansion rather than a spin-weighted spherical harmonic expansion.Using the inversion of the spin-weighted spherical harmonic expansion of the spinweighted spheroidal harmonics [58,59], one can re-expand Eq. ( 120) in terms of spin-weighted spheroidal harmonics, This completes our formalism as we have separated the source into the spin-weighted spheroidal harmonic modes of the Teukolsky equation.The radial Teukolsky equation (24) can now be solved with f ℓm [r] in (121).Our formalism is summarised in Fig. 3.

III. EXAMPLE: FIRST-ORDER TEUKOLSKY EQUATION WITH AN EXTENDED SOURCE
In this section, we present a toy-model application of our method to show that it provides a high-accuracy expression of the source of the Teukolsky equation.The accompanying Mathematica notebook [79] provides the implementation of this example for reference.This notebook was built upon other publicly available notebooks [49,80].
While the initial motivation for our method is separating the source of the second-order Teukolsky equation, it also applies to the more straightforward case of the first-order Teukolsky equation.For our toy model, we will apply our formalism to an extended source of the first-order Teukolsky equation Eq. ( 8).This is a sufficient example for applying our formalism because the second-order Teukolsky equation (Eqs.( 17) and (18), or the Campanelli-Lousto form [47]) are similar to the first order equation: the sources involve operators which depend on background Kerr quantities and the inputs can be expanded in terms of spin-weighted spherical harmonics.In upcoming papers, we will apply our formalism to the second-order self-force and quasi-normal mode problems.
Express the inputs A1, A2, ..., An in the form of Eq. ( 100), a spin-weighted spherical harmonic mode decomposition, up to factors of ρ and ρ.
Express f in GHP form using sections II B and II C. Apply the commutation relations, Eqs. ( 55) to (57) (and their complex conjugates).

Convert to Held form:
Replace the spin coefficients (and ψ2) with their Held form using Eqs.( 70) to (74); Convert the GHP derivatives to Held form using Eqs.( 66) to (68); Simplify the expression using Eqs.( 77) to (80); Convert to coordinate form: Replace Held scalars with their coordinate form using Eq. ( 102).
The resulting expression will be of the form Eq. (120).
Re-expand the source in terms of spin-weighted spheroidal harmonics using Refs.[58,59].The source is now in the form of Eq. ( 121) making the master Teukolsky equation separable, and one can now solve the radial Teukolsky equation, Eq. ( 24), for each ℓm mode.
FIG. 3. Summary of our method for separating the source of the Teukolsky equation.
For our example, we take Eq.(122) as our source and will express it in the form of Eq. (120).Following Fig. 3, we will expand the stress-energy tensor in terms of spinweighted spherical harmonics, similarly to Eq. (100).We also express Eq. ( 122) in the Held formalism as outlined in Sec.II D; please see the accompanying Mathematica notebook [79] for the step-by-step application and the resulting formulas.We require a stress-energy tensor with an extended radial profile.The stress-energy tensor of a perfect fluid is where ϱ is the density, p is the pressure, u a is the fourvelocity.To produce a simple example for applying our method, we assume p = 0 and where P and Λ are constants.Similarly, we assume the fluid is stationary in the rest-frame of the black hole, Using Eqs. ( 1), ( 59), ( 61), ( 123) and (125), where we have explicitly expressed the ρ and ρ content in Eqs. ( 126) to (128).Using Eqs. ( 103) and (104), we can decompose Eqs. ( 126) to (128) into spin-weighted spherical harmonics, In general applications of our formalism, the components of the stress-energy tensor will be expressed as a sum of spherical harmonic modes, in the form of Eq. (100).For our toy-model example, the simplicity of our stress-energy tensor has resulted in the mode decomposition being that of a single mode for each component (Eqs.(129) to (131)).Inserting Eqs.(129) to (131) into the Held formalism version of Eq. (122) allows us to extract all the ρ and ρ factors, expressing the source in the form of Eq. ( 105).The highest factor is ρ 6 ρ4 .We can pull out a common factor of ρ 6 ρ4 using Eq. ( 27) and combine the spin-weighted spherical harmonics using Eqs.( 35) and (37) to express Note there is no e −iωt dependence, and the spin-weighted spheroidal harmonics are equivalent to spin-weighted spherical harmonics because the source is stationary (ω = 0).
To complete our decomposition method requires expressing ρ 6 ρ4 in a separable form and combining it into the mode summation.To do so, we can use the Fourier expansion (with k = 8) for Σ −6 in Eq. ( 118), as ρ 6 ρ4 = Σ −6 (−r − ia cos[θ]) 2 .The resulting expression, after again applying Eqs. ( 35) and (37), is of the form The explicit form of fℓm [r] can be found in the accompanying Mathematica notebook.
In order to validate our expression for Eq. ( 132), we compare it to the four-dimensional expression for the source, Eq. ( 122).It is straightforward to calculate due to the simplicity of the stress-energy tensor using Eqs.( 54) and ( 58) to (61) (see Ref. [39] for expressions of the spin-coefficients in the Kinnersley tetrad).We compare our four-dimensional source to the sum in Eq. ( 132) in Fig. 4 near the horizon (r = r + ≈ 1.436M ).In the top plot, the two data sets overlap each other.The bottom plot shows that the error is less than 0.5% and takes the same form as the error in Fig. 1, as expected.One could decompose Eq. (134) into spin-weighted spherical harmonics using Eq. ( 25) (with spin-weighted spherical harmonics instead of spin-weighted spheroidal harmonics).Solving Eq. ( 25) requires computing an integral for each ℓm mode either analytically or numerically.For our toy-model example, the integrals are very slow to compute analytically; for second-order sources, which are much more complicated functions of r and θ (see Ref. [38] and the PerturbationEquations package in Ref. [59] for examples in Schwarzschild), analytical integration would be impractical.Calculating the integrals numerically requires computing an integral at each radial point on a grid, which is inefficient.Our formalism decomposes the source into modes without calculating integrals, so we expect it to provide efficiency savings for second-order calculations where the equations for the source are very long.
Near the horizon is the most inaccurate zone for our decomposition.This can be seen in Fig. 5, which plots the radial profile of the error for a = 0.9 and θ = 0.01.Additional precision errors can be encountered near the horizon if insufficient precision is used (less than 32 digits of precision).This is due to the Kinnersley tetrad being singular at the horizon.Hence, when expressing regular quantities using the Kinnersley tetrad near the horizon, cancellations of large numbers are required.This problem could be avoided by using a regular tetrad near the horizon, such as the Hartle-Hawking tetrad [69,73].Top: the source of the Teukolsky equation for a stationary ideal gas cloud is plotted in black (calculated using our decomposition method) and red (calculated explicitly) for a = 0.9, Λ = 1, P = 1, and r = r+ ≈ 1.436M .The lines overlay each other, showing our decomposition method is accurate.Bottom: the fractional disagreement between the two plots is shown; note the error is the same as Fig. 1, as expected.

IV. CONCLUSION
In this paper, we produced a formalism for decomposing the source of the Teukolsky equation into spin-weighted spheroidal harmonics analytically.Our formalism leverages the Held-formalism's [62] ability to extract the ρ and ρ dependency in background Kerr quantities.The remaining Kerr quantities exhibit a straightforward separation of variables in the Kinnersley tetrad [42].We then used a truncated Fourier series to separate the angular and radial dependency in ρ and ρ.An order k = 8 Fourier series has a less than 0.5% error for approximating Σ −6 .This error should be sufficiently small for second-order self-force calculations, as the error rapidly decreases as one moves away from the horizon, and second-order effects are small compared to first-order effects.Additionally, the error can be decreased by increasing the order of the The fractional error in the source of the Teukolsky equation for a stationary ideal gas cloud when calculated using our decomposition method.For a = 0.9, Λ = 1, P = 1, and θ = 0.01.The error is less than 0.5% and decreases as one moves away from the horizon as expected. expansion.
In section III, we demonstrated our method with a toy example: separating the source of the first-order Teukolsky equation for a stress-energy tensor of a pressure-less ideal gas.The calculation is publicly available in a Mathematica notebook [79], which explicitly goes through each step of applying our formalism, as summarised in Fig. 3.The calculations confirm the accuracy of our formalism, as we show the error in our source is the same as the error in our expansion for ρ and ρ, as expected.
There is an alternative method for separating the source of the second-order Teukolsky equation.Using the orthogonality of the spin-weighted spheroidal harmonics [50], as used in Ref. [29] for a quadratic quasi-normal mode calculation.However, our method provides the efficiency advantage of being analytic, avoiding the need to integrate at each radial point on a grid.This efficiency saving may be crucial for the more involved second-order self-force calculations in Kerr.
This paper focuses on applying our method in the frequency domain, but the formalism is also applicable in the time domain.In the time domain, decomposing into spin-weighted spherical harmonics reduces the 2+1 PDE Teukolsky equation into a coupled set of 1+1 PDEs with only nearest and next-to-nearest ℓ-mode coupling [86][87][88][89][90]. Hence, our formalism could help produce the first second-order self-force calculations in the time domain.

Take a tetrad related
to the Kinnersley tetrad where θ and A satisfy δθ = δA = 0; then β, β′ , and m a are only re-scaled by a factor of e iθ , so the relation of the GHP derivatives ð and ð ′ to the spin raising and lowering operators remains trivial.The Hartle-Hawking tetrad [69, 73] is one such example, where A = 2(r 2 +a 2 ) ∆ and θ = 0.However, this condition does not hold for all tetrads, such as the Carter tetrad.To transform from the Kinnersley to the Carter tetrad A = 2Σ ∆ and θ = −i ln ζ √ Σ .Hence,

FIG. 1 .
FIG.1.Top: the function Σ −6 plotted in blue (for a = 0.9 and r = r+ ≈ 1.436M ), and the Fourier expansion (for k = 8, see Eq. (114)) is over-plotted as a dashed yellow line.Bottom: the fractional disagreement between the two plots is shown.There is less than 0.5% error, and the error decreases for increasing radius (see Fig.2).

10 − 8
FIG.5.The fractional error in the source of the Teukolsky equation for a stationary ideal gas cloud when calculated using our decomposition method.For a = 0.9, Λ = 1, P = 1, and θ = 0.01.The error is less than 0.5% and decreases as one moves away from the horizon as expected.