Perturbations of spinning black holes in dynamical Chern-Simons gravity I. Slow rotation equations

The detection of gravitational waves resulting by the LIGO-Virgo-Kagra collaboration has inaugurated a new era in gravitational physics, providing an opportunity to test general relativity and its modifications in the strong gravity regime. One such test involves the study of the ringdown phase of gravitational waves from binary black-hole coalescence, which can be decomposed into a superposition of quasinormal modes. In general relativity, the spectra of quasinormal modes depend on the mass, spin, and charge of the final black hole, but they can be influenced by additional properties of the black hole, as well as corrections to general relativity. In this work, we employ the modified Teukolsky formalism developed in a previous study to investigate perturbations of slowly rotating black holes in a modified theory known as dynamical Chern-Simons gravity. Specifically, we derive the master equations for the $\Psi_0$ and $\Psi_4$ Weyl scalar perturbations that characterize the radiative part of gravitational perturbations, as well as for the scalar field perturbations. We employ metric reconstruction techniques to obtain explicit expressions for all relevant quantities. Finally, by leveraging the properties of spin-weighted spheroidal harmonics to eliminate the angular dependence from the evolution equations, we derive two, radial, second-order, ordinary differential equations for $\Psi_0$ and $\Psi_4$, respectively. These equations are coupled to another radial, second-order, ordinary differential equation for the scalar field perturbations. This work is the first attempt to derive a master equation for black holes in dynamical Chern-Simons gravity using curvature perturbations. The master equations can be numerically integrated to obtain the quasinormal mode spectrum of slowly rotating black holes in this theory, making progress in the study of ringdown in dynamical Chern-Simons gravity.


I. INTRODUCTION
The discovery of gravitational waves (GWs) has provided a new avenue for scrutinizing the predictions and phenomenology of Einstein's theory of general relativity (GR) in regimes characterized by non-linear and dynamic gravitational effects [1,2].GWs offer the opportunity to investigate the properties of astrophysical objects where gravity is notably intense, such as black holes (BHs) and neutron stars (NSs).In particular, GWs are often generated by the coalescence of binary BH systems, wherein two BHs orbit each other, gradually inspiraling due to the emission of GWs, and ultimately merging to produce a final BH that emits GW radiation as it settles down.GWs during this part of the coalescence, known as the ringdown phase, comprise a superposition of numerous quasinormal modes (QNMs), each with a complex eigenfrequency.By analyzing GW observations, it may thus become possible to efficiently explore the distinctive spectrum of QNMs exhibited by BHs as they ring down [3,4].
QNMs are the characteristic vibrational modes of BHs that are excited when the BH is perturbed.The study of QNMs can provide important information about the fundamental properties of BHs and their surrounding spacetime.In particular, the QNM spectrum of astrophysical (i.e., uncharged) BHs in GR is fully determined by just two parameters: the mass and spin of the remnant BH.One promising application of QNMs is to test modified gravity theories and alternative models of compact objects [5].The predictions of modified gravity theories may deviate from GR's and can manifest as modifications of the QNM spectrum of BHs [6][7][8][9][10][11][12].By observing the QNM spectrum of merging BHs with GW detectors, it may be possible to place constraints on these modifications to GR and its models given two or more detections of QNM frequencies with a strong signal-to-noise ratio.Recent GW detections of binary BH mergers have provided some of the most precise tests of GR in the strong-field regime [13].With improvements in detector technology and advancements in computational techniques, using ringdown to test modified gravity theories seems possible in the near future.
To study QNMs in GR, there are two well-established approaches within BH perturbation theory.One approach, proposed by Regge and Wheeler [14], perturbs the metric directly and it has been successfully applied to nonrotating and slowly rotating BHs in GR [14][15][16][17], as well as in modified theories of gravity [6][7][8][9][10][11][12].However, this approach has not been successfully applied to BHs with a general spin.An alternative approach uses curvature perturbations, and it was presented by Teukolsky [18] to study rotating Kerr BHs in GR, including their QNM spectrum and dynamical stability [19].This framework uses the Newman-Penrose (NP) formalism [20] and considers curvature perturbations characterized by quantities known as Weyl scalars.The success of the Teukolsky formalism lies in its ability to provide a decoupled evolution equation for each of the Ψ 0 and Ψ 4 Weyl scalars, which describe transverse gravitational perturbations, and physically represent ingoing and outgoing GWs, respectively.Not only are these quantities decoupled from other gravitational degrees of freedom, their evolution equations are also separable into a radial and an angular equation [4,19,21].
The Teukolsky formalism in GR requires the background geometry to be algebraically of type D under the Petrov classification [22], as is the case for Schwarzschild and Kerr spacetimes in GR.This Petrov type D property implies that four of the five Weyl scalars vanish on the background.However, when considering beyond-GR (bGR) theories, the deviations introduced may lead to BHs described by non-Petrov-type-D spacetimes.For instance, rotating BHs in dynamical Chern-Simons (dCS) gravity and Einstein-dilaton Gauss-Bonnet gravity are algebraically general and classified as Petrov type I.As a consequence, the Teukolsky formalism cannot be directly applied to bGR BH spacetimes.Therefore, calculating the QNM spectra of spinning BHs in bGR theories has been an open problem for a long time and warrants new approaches.One potential resolution to this problem became available recently with the development of the modified Teukolsky formalism [23][24][25][26].This formalism, in theory, enables studying perturbations of spinning BHs in bGR theories and calculating the QNM spectra of such non-Ricci-flat, matter vacuum Petrov type I BH spacetimes.Yet, for tests with GW data, a key theoretical challenge remains to calculate the QNM spectra of BHs in modified theories, and then compare them with observations.This work aims at computing the QNM spectrum of one type of BH spacetimes.
In this work, we restrict our attention to modifications to GR where a scalar field is non-minimally coupled to topological invariant quadratic terms in curvature.A subset of this class of theories, known as dCS gravity [27,28], was proposed to explain the matter-antimatter asymmetry of the universe.This is achieved by the introduction of additional parity-violating gravitational interactions, challenging a fundamental pillar of GR.Due to the quadratic nature of this theory, it is not strongly constrained using weak field tests and evades binary pulsar tests [29] and GW polarization tests [30].However, early work using the metric perturbation approach shows promise that ringdown tests can be useful in constraining this modified theory of gravity [9].
We focus on slowly rotating BHs in dCS gravity to leading order in spin in this work.This calculation serves as a validation of the newly developed formalism [23], as the results herein can then be confirmed with the results obtained using the metric perturbation approach [9,31].We first use the formalism prescribed in [23] to obtain the modified Teukolsky equation for a slowly rotating BH in dCS gravity.To leading order in spin, these BH backgrounds are described by a non-Ricci-flat, matter vacuum Petrov type D spacetime [23,32].Solving the Bianchi identities, we obtain the modified Teukolsky equation first in the null NP basis.We then rewrite this equation in the coordinate basis by defining a tetrad (similar to the Kinnersly tetrad in GR).Finally, we make use of the properties of spin-weighted spheroidal harmonics to eliminate the angular dependence and obtain a radial second-order differential equation.Due to the non-minimal coupling between the scalar field and the curvature, the perturbed master equations of the scalar field and the Weyl scalars (Ψ 0 or Ψ 4 ) are coupled in dCS gravity.Moreover, some of these quantities also require metric reconstruction within GR [33][34][35][36][37], making the problem more challenging, yet still solvable, as we demonstrate here.Having obtained the master equations in this work, we will use the eigenvalue perturbation method [24,38,39] in future work to calculate the QNM spectrum and compare it with the results using metric perturbations in dCS gravity [9,31].
The remainder of this paper is organized as follows.We first present in brief the action and the field equations of dCS gravity in Sec.II along with the slowly rotating BH solutions in this theory.In Sec.III, we present an overview of the modified Teukolsky formalism in [23], define a three-parameter expansion under the slow-rotation approximation, and calculate the NP quantities on the dCS BH background.In Sec.IV, we provide a concise review of the metric reconstruction procedures in GR.In Secs.V, VI, and VII, we calculate the source terms of the scalar field equation and the modified Teukolsky equation of Ψ 0 and Ψ 4 in the null NP basis, the results of which are summarized in Sec.VIII.In Sec.IX, we project the equations into the coordinate basis and extract their radial parts using the properties of spin-weighted spheroidal harmonics.Finally, in Sec.X, we summarize our work and discuss some future avenues.Henceforth, we adopt the following conventions unless stated otherwise: we work in 4-dimensions with metric signature (−, +, +, +) as in [40].For all NP quantities except the metric signature, we use the notation adapted by Chandrasekhar in [41].

II. BHS IN DCS GRAVITY
In this section, we will present the details of the theory and the background BH spacetime used in this work.

A. dCS gravity
In this subsection, we briefly review the dCS gravity following the discussion and the convention in [42].A more detailed review of dCS gravity can be found in [27,28].The action of dCS gravity is where κ g = 1 16πG , * R µνρσ is the dual of the Riemann tensor, * and ϑ is the pseudoscalar field coupled to the Pontryagin density P := R νµρσ * R µνρσ via the dCS coupling constant α.The quantity V (ϑ) is a potential for ϑ, which we set to zero for the reasons explained in [43][44][45], along with any matter contribution L matter (since we will work with matter vacuum BH spacetimes) for the remainder of this work.From Eq. ( 1), we find that [ϑ] = 1 and [α] = L 2 in geometric units.Using these coupling constants, we can define a dimensionless parameter ζ characterizing the strength of the dCS correction to GR, where ζ is defined in [42] to be with M the typical mass of the system.When the system under consideration contains a single black hole, then M is its mass.When considering a binary system, then different corrections to the solutions of the field equations will scale with different (dimension-4) combinations of the two masses.
We have here adopted the trace-reversed form of the field equations, using the fact that the C-tensor is traceless, as it will render future calculations simpler.The dCS action presented above is an effective theory that includes only linear in α and quadratic in curvature corrections to the Einstein-Hilbert action, thus ignoring higher order terms in α and in curvature.Therefore, the resulting field equations are also similarly effective, and their solutions ought to be truncated at leading order in α and considered only for systems (and regimes of spacetime) with small curvatures.Various previous work [42,46] have studied the regime of validity of this effective action and its curvature cutoff.In essence, the effective theory remains valid provided (α 2 /κ g )P 2 ≪ 1, where recall that P has been defined as the Pontryagin density.When this is the case, the higher order in α and in curvature terms neglected in the action above can continue to be ignored.The systems studied in this paper involve the exterior spacetime of remnant BHs with masses in the range 3M ⊙ < M < 10 7 M ⊙ .For such systems, the theory remains effective, and ζ ≪ 1 provided √ α ≪ 10 7 km, which will be assumed henceforth; the best current constraints on α come from NICER and advanced LIGO observations, and they require that α ≤ 8.5km at 90% confidence [47].In this range of α and for these BH masses, the quadratic curvature corrections to the Einstein Hilbert action will remain perturbative, and the higher-order in α and curvature terms will remain controlled relative to the quadratic term included in the dCS action.

B. Slowly rotating BHs in dCS gravity
Solutions to BHs in dCS gravity have been found both numerically [48] and analytically [42,[49][50][51].For this work, we consider the analytical solution found in [49,50], which were obtained by perturbatively solving the field equations ( 4) and ( 5) to linear order in both the dimensionless spin parameter χ, where χ ≡ a/M (with the dimensional spin parameter a = S/M , where S is the spin angular moment, and M is the BH mass), and the dCS expansion parameter ζ, defined in Eq. ( 3).This analytical solution casts the line element of a slowly rotating BH in dCS gravity as where, following the convention in [41], the line element for the Kerr metric in Boyer-Lindquist coordinates (t, r, θ, ϕ) is given by with ρ2 ≡ r 2 + a 2 cos 2 θ and ∆ ≡ r 2 − 2M r + a 2 .To leading order in χ and ζ, the dCS modification to the Kerr line element is given by where The above dCS correction to the line element is of O(ζ 1 , χ 1 , ϵ 0 ), and thus, we can use the tri-variate notation that we will introduce in Sec.III C to write The dCS metric of a slowly-rotating BH is then identical to the Kerr metric, except for the (t, ϕ) component, which acquires the correction presented above.Similarly, the background scalar field at O(ζ 1 , χ 1 , ϵ 0 ) is given by where we have absorbed a factor of ζ 1/2 into the expansion of ϑ as explained in more detail in [23].

III. BH PERTURBATIONS IN TEUKOLSKY FORMALISM
In this section, we review the modified Teukolsky formalism developed in [23].In this paper, we extend the twoparameter expansion scheme in [23] to a three-parameter expansion discussed in Sec.III C to incorporate the slowrotation approximation, following [9,31].

A. Modified Teukolsky Equation
As discussed previously in Sec.I, for studying perturbations of non-rotating BHs, we obtain the perturbed field equations and decompose these into master equations by making use of the metric perturbations [14][15][16]52].These metric perturbations are separated into two sectors, depending on their behavior under a parity transformation.For each parity, all the metric degrees of freedom are then packed into one master function: the Regge-Wheeler function for odd parity [14] and the Zerilli-Moncrief function for even parity [15,16].The master equations governing these master functions are decoupled from other dynamical degrees of freedom of the metric fields and are separable into radial and angular equations.
However, for rotating BHs in GR, due to the lack of spherical symmetry, to obtain the decoupled and separable perturbed field equations, one has to use the Teukolsky equations [18,19,53], where the fundamental variables to solve for are the Weyl scalars Ψ 0,4 characterizing curvature perturbations.In this case, the master equations of Ψ 0,4 are decoupled from other NP quantities and are separable into purely radial and purely angular equations.For a quick review of the NP formalism and the Teukolsky formalism in GR, one can refer to the original papers [18,19,53,54], the book [41], or more recent papers that work in the Teukolsky formalism [23,37].
In modified gravity, most calculations for rotating BHs have so far been done using metric perturbations and the slow-rotation expansion, e.g., [9,31] in dCS gravity, [10,11] in EdGB theory, and [12,55] in higher-derivative gravity.However, these approaches cannot deal with BHs with arbitrary spin, which motivated the development of the modified Teukolsky formalism in [23,24].Following the formalism in [23,24], one can find separable and decoupled equations for Ψ 0,4 of BHs with arbitrary spin in a wide class of modified gravity theories, such as in dCS gravity, which can be treated as an EFT extension of GR.In this paper, we will use the modified Teukolsky equations of Ψ 0,4 in [23].For an alternative approach following [56] by projecting the Einstein equations to the Teukolsky equations, one can refer to [24].
In [23], the authors introduced a two-parameter expansion, in terms of ζ and ϵ where 1. ζ is the parameter characterizing the strength of modifications to GR.In the case of dCS gravity, ζ is given by Eq. (3).
2. ϵ is the parameter characterizing the strength of gravitational perturbations, which also appears in GR.
In this way, one can expand all the NP quantities as and the extra non-metric fields, such as the pseudoscalar field ϑ in dCS gravity, as ϑ = ϑ (0) + ϵϑ (1) = ζϑ (1,0) + ζϵϑ (1,1) , where we have reserved the single superscript notation for only an expansion in ϵ.When using the double-superscript notation, however, the first superscript will also refer to contributions proportional to ζ to a given power.In contrast, the second superscript will refer to terms proportional to ϵ to a given power.Using the expansion in Eqs. ( 14) and ( 15), the authors in [23] found that for a rotating BH described by a matter vacuum, non-Ricci-flat, Petrov type I spacetime that perturbatively deviates from a Petrov type D spacetime in GR, the gravitational wave perturbation Ψ 0 satisfies where H (0,0) 0 is the Teukolsky operator for Ψ 0 in GR [18], and the source terms are divided into a "geometric piece", with

S
(1,1) and a "Ricci piece", with S 1,2 given by The operators H 0,1 , E 0,1 are defined as where Ψ 2 is a NP scalar, and we have also defined The operators that appear in the above definitions are defined for convenience to be where (D, ∆, δ, δ) are the usual NP differential operators (constructed by contracting the tetrad with partial derivatives), while (ε, ρ, µ, γ, α, β, π, τ ) are spin coefficients, with the overhead bar denoting complex conjugation, and (a, b, c, d) are certain constants.For a complete derivation of the equations above and the definition of Weyl scalars, spin coefficients, directional derivatives, and Ricci NP scalars, one can refer to [20,23].
For this work, we are only interested in studying slowly rotating BHs in dCS gravity up to O(χ).Such BHs are described by vacuum non-Ricci-flat Petrov type D spacetimes [32].The modified Teukolsky equations for these BHs hold the same form as given in Eq. ( 16) with the source terms S can be obtained from Eq. ( 16) by the Geroch-Held-Penrose (GHP) transformation [57] and is given in [23]: where is the Teukolsky operator in GR for Ψ 4 , and the "geometric piece" of the source terms is defined as with whereas the "Ricci piece" is defined as with The operators H 3,4 and E 3,4 are defined as with Although the formalism above works for BHs with arbitrary spin in dCS gravity, we choose to use the slowrotation expansion in this paper, so we can check the consistency of our results with prior work using metric perturbations [9,31] in our next paper [58].We implement a slow-rotation expansion of the above equations in Sec.III C.

B. Structure of the source terms
Here, we further discuss the structure of the source terms in Eq. ( 16) presented in [23].In particular, we will focus on the source terms that are non-vanishing for a non-Ricci-flat, Petrov type D BH in dCS gravity, given in Eqs.(18a) and (19).The source term in Eq. (18a) only depends on the perturbed Weyl scalar Ψ (0,1) 0 in GR and the dCS corrections to the stationary NP quantities at O(ζ 1 , ϵ 0 ).One can solve the Teukolsky equation in GR [18,19] to calculate Ψ (0,1) 0 directly.The NP quantities at O(ζ 1 , ϵ 0 ) can be computed from the dCS metric in Eq. (10), as shown in more detail in Sec.III D.
Due to the non-minimal coupling between the scalar field and the metric, the source term S (1,1) in Eq. (19) depends on both the scalar field perturbations and the metric perturbations.To compute S (1,1) , we first need to calculate the NP Ricci scalars Φ ij using the stress tensor or the Ricci tensor, i.e., where l µ is one of the NP tetrad basis vectors.Since the background and perturbed scalar field in GR vanish, we have that ϑ (0,0) = ϑ (0,1) = 0. Therefore, the NP Ricci scalars Φ (1,1) ij can be expressed as a function of the scalar field perturbation ϑ (1,1) and the metric perturbation h (0,1) as where g (0,0) and h (0,1) are shorthand for terms that depend on the metric tensor of the GR background and of the metric perturbation due to GWs reconstructed in GR, respectively.From Eqs. ( 19), (20), and (32), we notice that S (1,1) couples the GWs in GR and the scalar field ϑ, so we need to solve the equations of motions of these non-gravitational fields to find their contributions to the stress tensor and S (1,1) .Morevover, from Eqs. ( 21)-( 23), we see that S (1,1) in Eq. ( 19) depends on Ψ (0,1) 2 , the directional derivatives at O(ζ 0 , ϵ 1 ), and the perturbed spin coefficients at O(ζ 0 , ϵ 1 ), which need to be retrieved from the reconstructed metric perturbation h (0,1) µν for GR GWs.One can either follow the metric reconstruction approach in [33][34][35][36][59][60][61][62][63], the so-called Chrzanowski-Cohen-Kegeles (CCK) procedures, which involves defining an intermediate quantity called the Hertz potential, or the approach in [37,41], which solves the remaining NP equation directly.In this paper, we choose to follow the more widely used CCK procedures and apply them to compute the source term S (1,1) .

C. Slow-rotation expansion
When considering a slow-rotation expansion, in addition to the quantities given in Eqs. ( 14) and (15), one needs to consider an additional expansion in the dimensionless spin parameter χ = a/M .As an extension of Eqs. ( 14) and (15), all the NP quantities now admit a three-parameter expansion in ζ, ϵ, and χ, where as well as the pseudoscalar field, In what follows, it will sometimes be convenient to hide the χ expansion in more compact notation, such as ψ (1,1) = ψ (1,0,1) +χψ (1,1,1) .When only a two-parameter expansion is denoted, the χ expansion will be assumed.
The second set of terms are the GW perturbations to the Kerr metric in GR up to O(ζ 0 , χ 1 , ϵ 1 ).These terms include perturbed Weyl scalars, NP spin coefficients, and directional derivatives, all of which need to be evaluated in GR but include spin perturbations.To evaluate this type of terms, we need metric reconstruction of GW perturbations at O(ζ 0 , ϵ 1 ), the procedures of which are discussed in detail in Sec.IV.
The third set of terms are the one we want to solve for, which are corrections to GW perturbations in dCS gravity.The term Ψ (1,0,1) corresponds to gravitational perturbations sourced by non-rotating BHs in dCS gravity.Since ϑ (1,0,0) = 0, Ψ (1,0,1) is purely sourced by the leading contribution to the dynamical pseudoscalar field ϑ (1,0,1) , so only S (1,0,1) contributes to Eq. ( 16), and no metric reconstruction is needed.The term Ψ (1,1,1) corresponds to leading-order corrections to gravitational perturbations of slowly rotating BHs in dCS gravity.Unlike the nonrotating case, since both the metric and ϑ are corrected at can either be driven by dynamical GW perturbations in GR or dynamical ϑ.
For the first type of correction, the driving terms can come from S (1,1)   geo in the form of terms proportional to the product h (1,0) h (0,1) .As discussed in Sec.III A and [23], this kind of terms is due to the correction to the background geometry, so they are independent of bGR theories.Up to O(ζ 1 , χ 1 , ϵ 1 ), the background spacetime is still Petrov type D [42], so S 1,non-D = 0 in Eq. ( 17), and one does not need metric construction to evaluate these terms [23].In Sec.VI A, we will compute S (1,1)   geo in detail.Besides S (1,1)   geo , there is also contribution from S (1,1) in the form of terms proportional to the product ϑ (1,0) h (0,1) due to the effective stress tensor.In this case, metric reconstruction is needed, and we will compute S (1,1) in Sec.VI B.
In Sec.VI B, we will compute the source terms driven by ϑ (1,1) but leave ϑ (1,1) unevaluated.In our follow-up work [58], we will solve both the modified Teukolsky equation and the scalar field equation jointly to find the QNM shifts.Since metric reconstruction at O(ζ 0 , ϵ 1 ) is required for both the modified Teukolsky equation and the scalar field equation, we present a review of the procedures in Sec.IV.

D. NP quantities on background
In this subsection, we will present the background tetrad for a non-Ricci-flat, Petrov type D, slowly rotating, dCS gravity BH spacetime.To obtain the null tetrad for the metric given by the line element in Eq. ( 8), one can follow the general procedures prescribed in [23] or the standard procedures for finding the Kinnersley tetrad in GR given in [41].In [64], such a tetrad was found by following the second approach.Nonetheless, for completeness, let us re-derive the tetrad following the prescription in [41,64].Our result is consistent with the one in [64], but with additional tetrad rotations to set Ψ (1,0) 0,1,3,4 = 0. To begin, we first find the null geodesics in the equatorial plane for a dCS BH to be where E = −∂L /∂t is a constant of motion, L is the Lagrangian for Kerr in [41], ∆(r) = ∆(r) + 2aM G(r)/r + G(r) 2 /4, and G(r) = ζχ G(r).Following the procedures outlined in [41], we align the tetrad basis vectors l µ and n µ along the outgoing and ingoing null geodesics respectively at the equilateral plane with E = 1 such that where N is the normalization factor introduced to impose l µ n µ = −1.Since l µ and n µ are along null geodesics, l µ l µ = n µ n µ = 0 is satisfied automatically.Expanding Eqs.(38) and (39) up to O(ζ 1 , χ 1 , ϵ 0 ), we find where r s is the Schwarzschild radius given by r s = 2M , and Ñ (r) = (r − r s )/2r.When ζ = 0, l µ and n µ reduce to the Kinnersley tetrad of Kerr BHs expanded to O(χ 1 ).The tetrad basis vectors l µ and n µ in Eqs.(40) and (41) are the same as the principal null directions in Eq. ( 31) of [64].
To obtain the remaining components of the null tetrad, notice that the correction to the Kerr metric due to dCS gravity enters at O(ζ 1 , χ 1 , ϵ 0 ) only in the tϕ-component.Therefore, it can be expected that the corrections to the Kinnersley tetrad are only along the ∂ t and ∂ ϕ directions, which is seen to be true for l µ and n µ .Thus, at O(ζ 1 , χ 1 , ϵ 0 ), the corrections to the remaining null tetrad components, m µ and mµ , take the form where mµ can be obtained by taking the complex conjugation of Eq. ( 42).Imposing the remaining orthogonality Notice that m µ (and therefore mµ ) holds the same form as the Kinnersley tetrad of Kerr BH expanded to O(χ).

IV. METRIC RECONSTRUCTION
This section reviews how to reconstruct the perturbed metric and the corresponding NP quantities from solutions to the Teukolsky equation for Kerr BHs in GR.There are two approaches to metric reconstruction in general: the first approach involves systematically solving the Bianchi identities, Ricci identities, and commutation relations [37,41], whereas the second approach, or the CCK procedures, utilizes an intermediate Hertz potential to reconstruct the metric [33][34][35][36][59][60][61][62][63].In this work, the second approach is employed to perform metric reconstruction.

A. Metric perturbations
In this subsection, we present the reconstructed metric perturbation h (0,1) µν for GR GWs.For convenience, in this section, we will drop the superscript (0, 1) of h (0,1) µν and always assume that h µν is at O(ζ 0 , ϵ 1 ).The CCK procedures can be carried out in two different gauge choices: Ingoing radiation gauge (IRG): Outgoing radiation gauge (ORG): where h is the trace of h αβ with respect to the background metric.The reconstructed metric h αβ in the IRG and ORG are given in Eqs. ( 48) and ( 49), respectively [35,36], where the notation for the derivatives is given by Eq. ( 23), and ΨH is the complex conjugate of the Hertz potential.
We have also dropped the superscript (0, 1) of the Hertz potential for simplicity.Notice, since we use an opposite signature from [35,36], our Eqs.( 48) and ( 49) have an opposite sign from the results in [35,36].
Substituting the differential operators D † mω and D mω in Eq. ( 55) into the expression for the radial part of the Hertz potential in Eq. ( 54), we have where we have made use of the radial Teukolsky equation to reduce all second-and higher-order derivatives of s R (0,1) ℓm (r).In Eq. ( 57), the prime denotes the first derivative with respect to the radial coordinate r.The functions f i are spin weight s and mode dependent.These functions are lengthy and non-illuminating, so they have been presented in a separate Mathematica notebook [68].

C. Spin-weighted spheroidal harmonics
Spin-weighted spheroidal harmonics that appear in Eq. ( 53) are solutions to the angular Teukolsky equation in GR [18,19].In general, these are eigenfunctions of an equation of the form [69] where s represents the spin weight, s A γ ℓm is the eigenvalue of the equation, which has been numerically calculated in the literature [19].Comparing Eq. ( 58) with the angular Teukolsky equation in GR [18], we see that γ = χM ω.In the slow-rotation expansion, spin-weighted spheroidal harmonics s Y γ ℓm can be expanded as [19,69] where s Y ℓm are spin-weighted spherical harmonics with spin weight s.The factors s b m ℓ,ℓ±1 in Eq. ( 59) hold the form [69] To evaluate spin-weighted spherical harmonics, one can use In the special case that s = 0, the spin-weighted spherical harmonics become the standard spherical harmonics, i.e., Spin-weighted spherical harmonics and spin-weighted spheroidal harmonics also obey the following orthogonality relations, where dS is the solid angle element, and the integration is over the entire 2-sphere.At certain places, we might drop the superscript γ of s Y γ ℓm denoting its eigenvalue for simplicity.

D. Perturbed NP quantities
As we are working within the NP basis, in addition to the perturbed metric given by Eq. ( 48), we also require the reconstruction of the perturbed NP quantities, such as the perturbed tetrad, Weyl scalars, and spin coefficients.We adopt the methodology outlined in [37,65] to perform this reconstruction.The first step involves expressing the perturbed tetrad in terms of the background tetrad.This is accomplished by expanding the perturbed tetrad in terms of the background tetrad, and then utilizing the transformation properties of the tetrad to obtain the perturbed tetrad components in terms of the background tetrad components such that where e µ a represents a null tetrad such that and A a b are coefficients that map the background tetrad to the perturbed tetrad.
As shown in [37,65], one can always perform tetrad rotations to set six real parameters of the A a b(0,1) coefficients to zero.Then, expanding h µν in terms of e a(0,1) µ and e a(0,0) µ and using the linearized completeness relation, we find Comparing Eq. ( 70) to Eq. ( 68), we find [37,65], where we have dropped the superscripts of e µ(0,0) a for simplicity.Since we have adopted the sign convention in [41], our signature is opposite to that used in [37,65].Therefore, the perturbed tetrad in Eq. ( 71) has an opposite sign from the results of [37,65], as seen in Eq. (70).Equation ( 71) works for both the IRG and ORG.In the IRG or ORG, we can further set h m m = 0 in Eq. (B4), respectively, where a is any index in the NP basis.
For the spin coefficients, one can linearize the commutation relation following [41], where γ a bc is the Ricci rotation coefficients.Using the relation between spin coefficients and Ricci rotation coefficients in Eq. (B1), one can write C ab c in terms of spin coefficients, as listed in Eq. (B2) [20,41].From Eq. (B2), one can also solve for spin coefficients in terms of C ab c , and the results are in Eq. (B3).Expanding Eq. (72) using Eq. ( 68), one finds where the superscripts of C ab c(0,0) are dropped for convenience.The coefficients A a b(0,1) can be retrieved from Eq. ( 71).The GR structure constants C ab c(0,0) are directly given by Eq. (B2) and the spin coefficients in GR.With all the quantities in Eq. ( 73), one can then use Eq. ( 73) and (B3) to evaluate the spin coefficients at O(ζ 0 , ϵ 1 ).We have listed our result in Eq. (B4), which works for both the ORG and IRG.Our result is consistent with [37] up to the overall minus sign due to different signatures, which corrects some errors in [65].
To reconstruct Weyl scalars, one can either directly linearize the Riemann tensor and project it onto the NP basis to find Weyl scalars or use the Ricci identities in Eq. (B5).For both approaches, we use the perturbed tetrad in Eq. ( 71), and we check that the results are consistent.We also compare our results in the IRG to the equations in [35], which corrected a factor of 1/2 missed in [59] and are listed in Eq. (B6).After expressing everything in terms of the Hertz potential, our results of Ψ (0,1) 0,1,2,4 in the IRG agree perfectly with Eq. (B6) but not for Ψ (0,1) 3 .Since Ψ (0,1) 3 is not invariant either under tetrad rotations or coordinate transformations at O(ζ 0 , ϵ 1 ), this difference indicates that we might have a O(ζ 0 , ϵ 1 ) difference in the choices of coordinate or tetrad.
V. THE EVOLUTION EQUATION FOR ϑ (1,1) IN THE IRG In this section, we project the equation governing ϑ (1,1)  [Eq.(36)] onto the NP basis using the IRG.For convenience, we define the right-hand side of Eq. (36) as so that the evolution equation for ϑ (1,1) becomes This equation is first expressed in terms of the NP quantities, following which we evaluate its left-hand side using the background NP quantities in Sec.III D and Appendix A and its right-hand side using the reconstructed NP quantities at O(ζ 0 , ϵ 1 ) in Sec.IV and Appendix B.
The same methodology demonstrated in this section is applied to computing the modified Teukolsky equation in Sec.VI. Figure 1 presents a schematic illustration of the steps involved in calculating a completely separated radial evolution equation for the scalar field perturbation in the IRG.

A. Projection onto the NP basis
In this subsection, we project Eq. ( 36) onto the NP basis.This projection involves the projection of two fundamental quantities: the D'Alembertian operator □ and the Pontryagin density R * R onto the NP basis.In particular, our goal is to express these quantities in terms of NP quantities, particularly the below-mentioned quantities in the NP basis, namely Here, η ab is the metric in the NP basis.The notation f ,a denotes the directional derivative of f defined by the tetrad basis e µ a .The quantities γ abc are Ricci rotation coefficients, which can be expressed in terms of spin coefficients using Eq.(C1) [20].The tensor R abcd can be expressed in terms of Weyl scalars using Eq.(C3) [20].Therefore, the Pontryagin density and the □ operator can be expressed in the NP basis as The factor of i in Eq. (79a) arises from the normalization of the Levi-Civita tensor ϵ abcd in the NP basis.In the literature, such as in [71], the covariant Levi-Civita tensor is typically defined as denotes the Levi-Civita symbol.However, this definition encounters issues when attempting to convert a Levi-Civita tensor from Boyer-Lindquist coordinates to the NP basis, due to the determinant of the Jacobian relating these two bases often being complex.Thus, to convert the tensor density εµν•••γ to a tensor, we instead need to define which has the same normalization factor as the Einstein-Hilbert action.The absolute value in the usual definition is to impose that the Levi-Civita tensor is a real tensor in the Lorentzian signature, but the minus sign of √ −g will do the same trick.Since η = 1, we find the normalization factor to be i rather than 1 from Eq. ( 80).This is also consistent with that which shows that ϵ lnm m is an imaginary number.We have now expressed all the terms in Eq. ( 36) in the NP basis.
B. Left-hand side of Eq. (36) In this subsection, we compute the operator □ (0,0) acting on ϑ (1,1) to obtain the homogeneous component of Eq. (36).The operator □ (0,0) can be evaluated directly using the Kerr metric presented in Eq. ( 9).Alternatively, one can use Eq.(79b) and the NP quantities of Kerr, expanded up to O(χ), as given by Eqs.(A4)-(A6), and then setting ζ to zero.We therefore find where H (0,0) ϑ is the Teukolsky operator for particles with spin weight s = 0 in [18], where we have only kept the terms up to O(χ) and separated ϑ (1,1) as or in the slow-rotation expansion (85) Thus, in the absence of sources, Θ ℓm (r) satisfies Scalar field equation for ϑ (1,1) in the IRG Left-hand side of Eq. ( 36) Right-hand side of Eq. ( 36) In coordinate basis Equation (83) Separated radial equation Equation ( 86) Source term S (1,1) ϑ rewritten as Eq. ( 91) In coordinate basis using the Hertz potential Eq. ( 92) Expression for (R * R) (0,1) and □ (0,0,1) ϑ (1,1,0) in {t, r, θ, ϕ} coordinates in Eq. ( 93) and Eq. ( 95) Using the recipe in Sec.IX A Source terms as functions of the radial coordinate in the IRG, V R ℓm (r) and V 2 ℓm (r), are given by Eqs. ( 133) and ( 134) The radial evolution equation of the scalar field perturbation ϑ (1,1) for slowly rotating BHs in dCS gravity in the IRG is given by Eq. ( 138) FIG. 1.A schematic flowchart prescribing the steps involved in obtaining a separated radial evolution equation for the scalar field perturbation ϑ (1,1) for slowly rotating BHs in dCS gravity.This flowchart summarizes the details of the calculations described in Sec.V and Sec.IX B. Initial and final outcomes are represented by rectangular boxes, while intermediate results are symbolized by encapsulating bubbles.The directional arrows are meant to seamlessly guide the reader through the logical flow of the calculations.
where 0 A ℓm is the Teukolsky's separation constant for s = 0 [18].We therefore see that the left-hand side of Eq. ( 36) is separable in the radial and angular coordinates.Further, in Sec.IX, we show that the complete expression in Eq. ( 36) can be separated into radial and angular parts using spin-weighted spheroidal harmonics of spin weight zero.
For the first piece, according to Eq. ( 87), we essentially need to evaluate Ψ  53).Since we use the slow-rotation approximation in this work, we can further decompose spin-weighted spheroidal harmonics in terms of spin-weighted spherical harmonics using Eqs.( 59) and ( 60) such that Eq. ( 53) becomes Now, one can insert into Eqs.(B6c) and (87) the decomposition in Eq. ( 92) and the background NP quantities at O(ζ 0 , ϵ 0 ) in Eqs.(A4)-(A6) after setting ζ = 0.After using Eqs.(67a) and (67b) to simplify the terms with δ (0,0) and δ(0,0) acting on s Y ℓm (θ, ϕ), we find where f ′ (r) denotes the derivative of f with respect to the r coordinate for any function f (r).Here, 2 Rℓm (r) is the radial function of the Hertz potential for slowly rotating Kerr BHs in GR, which can be computed from the radial function of Ψ (0,1) 0 using Eq.(54a).One can, in principle, expand 2 Rℓm (r) further in χ and drop additional terms above O(χ) in Eq. (93).For simplicity, we choose not to do this additional expansion here but implement it when computing QNMs in [58], where we need to explicitly evaluate the radial functions of Ψ (0,1) 0,4 .We have also used the radial Teukolsky equation to reduce any n-th radial derivative of 2 Rℓm (r) with n ≥ 1 to 2 Rℓm and 2 R′ ℓm (r).The explicit forms of g ℓm i (r, ω, M ) are long and therefore presented through a separate Mathematica notebook as Supplementary Material [68].
For the second piece, according to Eq. ( 91), there is only a contribution from □ (0,0,1) ϑ (1,1,0) since the scalar field at O(ζ 1 , χ 0 , ϵ 0 ) vanishes in dCS gravity.Thus, we only need the reconstructed metric at (ζ 0 , χ 0 , ϵ 1 ).Using Eqs.(92), (A6), and (48), we find h (0,0,1) Now, we can evaluate all the directional derivatives and spin coefficients at O(ζ 0 , χ 0 , ϵ 1 ) using Eqs.(B4), (B8), and (B10).In the end, using Eq. ( 90), we find where ϑ R (r) is the radial part of the background scalar field in Eq. ( 13).In Eq. ( 95), the reconstructed metric only has contribution at O(χ 0 ), so the radial function 2 Rℓm (r) is evaluated on the Schwarzschild background.Since we choose not to expand R(r) in χ here, we do not distinguish 2 Rℓm (r) evaluated on the Schwarzschild or slowly rotating Kerr background.Combining Eqs. ( 93) and (95), we have the source terms in the equation of ϑ (1,1) up to O(χ).The master equation of ϑ (1,1) in the IRG in the coordinate {t, r, θ, ϕ} is presented in Eq. (125) of Sec.VIII.In Sec.IX, we will show that Eq. ( 76) up to O(χ) can be separated into a radial and an angular equation.In the following section, we apply the same procedures to evaluate the source terms in the modified Teukolsky equation of Ψ (1,1) 0 in terms of the reconstructed NP quantities and project the equation into the coordinate basis.

IN THE IRG
In this section, we evaluate the modified Teukolsky equation of Ψ (1,1) 0 in Eq. ( 16) following the similar pro-cedures in Sec.V. We first evaluate the left-hand side of Eq. ( 16) and the source term S (1,1)   geo due to the correction to the background geometry using the background NP quantities in Sec.III D and Appendix A. We then project the source term S (1,1) onto the NP basis and compute its coordinate-based value using the reconstructed NP quantities in the IRG given in Sec.IV and Appendix B. Figure 2 presents a schematic illustration of the steps involved in calculating a completely separated radial evolution equation for the Ψ (1,1) 0

Weyl scalar perturbation in the IRG.
A. Right-hand side of Eq. ( 16) and S (1,1)   geo Since slowly rotating BHs in dCS gravity are Petrov type D up to O(ζ 1 , χ 1 , ϵ 0 ), we only need to compute the Teukolsky operator H (0,0) 0 in GR and its stationary correction H (1,0) 0 in Eq. (18a).Thus, we do not need metric reconstruction in this subsection.
Using the Weyl scalars and spin coefficients found in Sec.III D and Eq. ( 21), we find that Perturbed gravitational field equation for Ψ (1,1) 0 in the IRG Left-hand side of Eq. ( 16) Right-hand side of Eq. ( 16) In coordinate basis Equation (96)

Separated radial equation
Left-hand side of Eq. ( 146) Source term S (1,1) Source term S (1,1)   geo In coordinate basis using the Hertz potential Eq. ( 92) terms proportional to Θ ℓm (r) S The radial evolution equation of the Weyl scalar perturbation Ψ (1,1) 0 for slowly rotating BHs in dCS gravity in the IRG is given by Eq. ( 146) FIG. 2. A schematic flowchart prescribing the steps involved in obtaining a separated radial evolution equation for the gravitational field perturbation described by the Ψ (1,1) 0 Weyl scalar in the IRG for slowly rotating BHs in dCS gravity.This flowchart summarizes the details of the calculations described in detail in Sec.VI and Sec.IX C. Initial and final outcomes are represented by rectangular boxes, while intermediate results are symbolized by encapsulating bubbles.The directional arrows are meant to seamlessly guide the reader through the logical flow of the calculations.
B. Source term S (1,1)   In this subsection, we evaluate the source term S (1,1)  from the effective stress tensor or the source term of the trace-reversed Einstein equation in Eq. ( 4).The first step is to project the Ricci tensor R µν onto the NP basis such that we can express the NP Ricci scalars Φ ij , where in terms of NP quantities (Weyl scalars, spin coefficients, and tetrad) and the scalar field ϑ.The precise relation between Φ ij and R ab is given in Eq. (C6).Using Eq. ( 4), we find where we have inserted an additional ζ − 1 2 to the term linear in ϑ and an additional ζ −1 to the term quadratic in ϑ to compensate the factor of ζ 1 2 we have absorbed into the expansion of ϑ in Eq. (15).Since ϑ enters at least at O(ζ), all the metric fields at the right-hand side of Eq. ( 103) are at O(ζ 0 ), which can be expressed in terms of NP quantities in GR.
Since we are interested in gravitational perturbations of vacuum spacetime, R (0,0) µν = R (0,1) µν = 0, and all the metric fields in Eq. (103) enter at O(ζ 0 ), the first term in Eq. ( 103) vanishes.Evaluating the rest of the terms, we find the seven independent components of R ab in terms of Weyl scalars, spin coefficients, directional derivatives, and the scalar field ϑ in Eq. (C4).
We can now evaluate the source terms in the modified Teukolsky equations.Inspecting the source term S (1,1) in Eq. ( 19), we can divide it into two parts based on whether S 1,2 are dynamical, For S (1,1) I , one can directly evaluate S (1,0) 1,2 in terms of the stationary scalar field ϑ (1,0) and the metric in GR using Eqs.( 20), (C4), and (C6).Then we only need to reconstruct the operators E 1,2 , we can further divide them into two parts based on whether Φ ij are dynamical, where we used that κ (0,0) = σ (0,0) = λ (0,0) = 0. Based on this classification, we can then additionally separate S

IIB , S
(1,1) For S (1,1) IIA , similar to S (1,1) I , we only need to evaluate Φ (1,0) ij using the background metric and the stationary scalar field ϑ (1,0) .The only quantities need metric reconstruction are these additional O(ζ 0 , ϵ 1 ) operators acting on Φ 02 .To organize our calculations, we can divide Φ (1,1) ij into four parts based on which kind of terms need metric reconstruction, similar to what we have done in Sec.V. Inspecting Eq. (C4), we notice that all the terms are some coupling of a Weyl scalar, a scalar field, a spin coefficient, and directional derivatives, which is due to the structure of R µν in Eq. ( 103).In this case, we make the following classification: 1. Weyl scalars are at O(ζ 0 , ϵ 1 ).In this case, the scalar field ϑ is at O(ζ 1 , ϵ 0 ).As discussed in Sec.III C, the leading contribution to ϑ (1,0) is at O(ζ 1 , χ 1 , ϵ 0 ) since non-rotating BHs in dCS gravity are still Schwarzschild.Then since we are interested in O(ζ 1 , χ 1 , ϵ 1 ) corrections, all the spin coefficients and directional derivatives are at O(ζ 0 , χ 0 , ϵ 0 ), the order of the Schwarzschild background.
3. Tetrad/directional derivatives are at O(ζ 0 , ϵ 1 ).Similar to the first two cases, since ϑ is at O(ζ 1 , χ 1 , ϵ 0 ), so all the Weyl scalars and spin coefficients can be evaluated on the Schwarzschild background.
4. The scalar field ϑ is at O(ζ 1 , ϵ 1 ), which has contributions from both O(ζ 1 , χ 0 , ϵ 1 ) and O(ζ 1 , χ 1 , ϵ 1 ) terms.Then, all the NP quantities generally need to be evaluated on the Kerr background expanded in the slow-rotation expansion to O(χ 1 ).Since ϑ (1,1)  also requires us to solve the scalar field equation in Eqs.(36), we choose not to compute ϑ (1,1) in this work but only list the source terms in terms of it.
We will solve the scalar field equation together with the modified Teukolsky equation in our follow-up work [58].
At O(ζ 1 , χ 1 , ϵ 1 ), using the classification above, we can set many terms to 0 since they are evaluated on the Schwarzschild background (i.e., when ϑ is stationary).Similar to Sec.V, the results of Φ up to O(ζ 1 , χ 1 , ϵ 1 ) are expressed in terms of the perturbed Weyl scalars, metric perturbations, and dynamical ϑ in Appendix C. Due to the complication of S (1,1) , we will not present the results here but provide its coordinate-based values directly in the next subsection.
For the first part, we find its coordinate-based value to be S S(1,1) In Sec.IX C, after getting the radial part of the equation of ϑ (1,1) , we will further express Θ ′′ R (r) in terms of Θ ℓm (r), Θ ′ ℓm (r), 2 Rℓm (r), and 2 R′ ℓm (r).Similarly, we find the second part to take the form S(1,1) where 2 Rℓm (r) is the radial function of the Hertz potential given by Eq. (54a), and 2 Rℓm (r) is the complex conjugate of 2 Rℓm (r).The radial functions q ℓm i (r, ω, M ) and qℓm i (r, ω, M ) are presented in a Mathematica notebook as Supplementary Material [68].Note that we have used the radial Teukolsky equation to eliminate any beyondfirst-order derivatives of the Hertz potential to obtain a simplified expression in Eqs.(109a) and (109b).The master equation of Ψ (1,1) 0 in the IRG in the coordinate {t, r, θ, ϕ} is presented in Eq. ( 125) of Sec.VIII.

IN THE ORG
In this section, we evaluate the equations governing the evolution of the perturbed scalar field ϑ (1,1) and the perturbed Weyl scalar Ψ in the IRG, as briefly discussed in Sec.IV B, the evaluation is more convenient in the ORG.We follow a set of steps similar to those in the IRG for ϑ (1,1) in Sec.V and for Ψ (1,1) 0 in Sec.VI.Below, we present the master equations of ϑ (1,1)  and Ψ (1,1) 4 in the ORG.
A. The equation of ϑ (1,1)   The scalar field perturbations are governed by Eq. (36).We now represent the right-hand side of Eq. ( 36) as Projecting the Pontryagin density onto the NP basis leads to the same set of equations as described in Sec.V A since our choice of gauge does not affect the quantities shown in Eqs.(79).
The master equation for the scalar field perturbations in the ORG are same as that shown in the IRG with H (0,0) ϑ and ϑ (1,1) both given in Eqs. ( 83) and ( 84) respectively, whereas T (1,1) ϑ is given in Eq.( 110).The lefthand side of Eq. ( 111) in the ORG remains unchanged from the IRG since the operator acting on the scalar field perturbations is evaluated on the background.On the other hand, since the source term in Eq. ( 111) depends on perturbed quantities, the value of these quantities is gauge dependent.In the ORG, the Pontryagin density given in Eq. (87) holds the following form in the coordinate basis where functions g i (r, ω, M ) are presented in a separate Mathematica notebook as Supplementary Material [68], and is the radial function of the Hertz potential for a slowly rotating Kerr BHs in GR computed from the radial function of the Ψ (0,1) 4 using Eq.(54b).To evaluate the remaining part of the source term T (1,1) ϑ , we use the perturbed spin coefficients given in Eq. (B3) and metric perturbations given in Eq. (49).We obtain where the functions h i (r) are presented in a separate Mathematica notebook as Supplementary Material [68], and ϑ R (r) is the radial part of the background scalar field given in Eq. ( 13).Similar to the case evaluated in the IRG in Sec.V, the radial function −2 Rℓm (r) is evaluated on the Schwarzschild background.Combining Eqs. ( 113) and ( 112) gives us the complete source term T (1,1) ϑ . The master equation of ϑ (1,1) in the ORG in the coordinate {t, r, θ, ϕ} is presented in Eq. (125) of Sec.VIII.

B. The equation of Ψ (1,1) 4
In this subsection, we present the modified Teukolsky equation for the Weyl scalar perturbation Ψ (1,1) 4 given in Eq. ( 24) in the coordinate basis.Following steps similar to Sec.VI, we separate the source terms into two categories: T (1,1)   geo and T (1,1) whose forms in the NP basis have been given in Eqs. ( 25)-( 27).
(118) For the same reason in Sec.VI A, we do not need metric reconstruction to compute T (1,1)   geo since the background spacetime is still Petrov type D. The source terms T (1,1)   geo hold the form where we have absorbed the factor of ρ 4 into H (1,1,0) 4 , and the functions D i (r) are presented in Appendix A.
2. T (1,1)   Using the expression for the metric perturbation in the ORG given in Eq. ( 49), one can evaluate the perturbed spin coefficients and perturbed Weyl scalars at O(ζ 0 , ϵ 1 ).These can then be used to evaluate the following source terms, which can be divided into two parts based on whether S 3,4 are dynamical.
Analogous to Sec.VI, T consists of terms dependent on the stationary scalar field, the background, and the perturbed metric in GR.
Similarly, we can further divide T

VIII. EXECUTIVE SUMMARY OF ALL MASTER EQUATIONS
This section presents an executive summary of the main results of this paper, whose derivation was presented in Secs.V, VI and VII.Through Secs.V, VI and VII, we used the tetrad defined in Sec.III D to rewrite Eqs. ( 16) and (36) using the IRG and Eqs. ( 24) and (36) using the ORG.In this section, we summarize these results and condense them into a single master equation for convenience.In later sections, we will present and apply a procedure to decouple the master equations for the scalar field perturbation ϑ (1,1) and the Weyl scalar perturbations Ψ (1,1) 0,4 ) into a set of coupled radial ordinary differential equations in the IRG (or ORG).
With this in mind, the master equations for ϑ (1,1)  [Eq.(36)], Ψ (1,1) 0 [Eq.( 16)], and Ψ (1,1) 4 [Eq. ( 24)] can all be expressed as where H represents the Teukolsky operator for a spin s field in GR given in [18].Recall that χ is the dimensionless spin parameter, and M is the BH mass.In Table I, we present the field quantities ψ which satisfy this equation and the source terms, ξ(r) and S on the right-hand side of Eq. ( 125), dependent on the gauge and the spin weight s of these fields.Observe that, clearly, the differential operator H in Eq. ( 125) acting on the field quantity ψ is exactly the same as the one derived by Teukolsky in [18] in GR for Kerr BH perturbations [c.f.Eq. (4.7) therein], expanded to leading order in spin.In addition, from Table I and the source terms in Eqs. ( 101), ( 108), ( 109), ( 119), (123), and (124), we notice that the (l, m) and (l, m) modes of Ψ (1,1) 0,4 need to be solved jointly, as we will discuss in more detail in Sec.IX.

IX. SEPARATION OF VARIABLES AND EXTRACTION OF THE RADIAL MASTER EQUATION
In this section, we extract the radial parts of the master equations for ϑ (1,1) [Eq.( 36)] (both in the IRG and ORG), for Ψ (1,1) 0 [Eq.( 16)], and for Ψ (1,1) 4 [Eq. (24)].Let us first present our procedures for eliminating the angular dependence in these equations and then apply them to specific cases.

Equation (
On the other hand, in GR, one has the well-known symmetry [21] so both Eqs.(126a) and (126b) contain source terms proportional to e −iω ℓm t and e −iω ℓ−m t .Thus, one has to consider Eqs.(126a) and (126b) jointly or solve the linear combination , i.e., Ψ (0,1) Ingoing radiation gauge Outgoing radiation gauge s (spin weight) TABLE I.In this table, we present the quantities Ψ, the spin weight s and the source terms S that appear in Eq. ( 125).
where we have absorbed an overall factor into the normalization of s R (0,1) ℓm (r) and s ψ ℓm (r).Furthermore, we have also taken ω ℓ−m not just in GR but also in dCS gravity due to the structure of Eq (126).As discussed and shown in more detail in [72], one can solve for both the complex constant η ℓm and the QNM frequencies ω ℓm using the eigenvalue perturbation method method in [24,38,39].The combination (η ℓm , ω ℓm ) have two independent solutions [24,25,72], resulting in the breaking of isospectrality.In this case, by plugging the ansatz in Eq. ( 128) into the scalar field equation or the modified Teukolsky equations, using Eq. ( 126), and matching the phase of the terms, we find where we have divided a factor of η ℓm in Eq. (129b), and s ψ ℓ±m (r) are radial solutions tied to (η ℓm , ω ℓm ).Similar procedures for generic modified gravity theories can be found in [72].
Notice that s f k ℓm (θ) are angular functions in the modified Teukolsky equation for the particle of spin weight s and mode (ℓ, m).The subscript s and subscripts (ℓ, m) do not indicate the mode number of the angular function itself.For example, 0 f k ℓm (θ) contains terms proportional to sin θ ±1 Y ℓm (θ, ϕ).
6. Since we use the slow-rotation approximation in this work, when computing the integrals involving s Y ℓm (θ, ϕ), one can further expand s Y ℓm (θ, ϕ) in terms of s Y ℓm (θ, ϕ) using Eq. ( 59).Thus, there are only spin-weighted spherical harmonics in these integrals.
7. After the angular integration, the angular functions s f ℓm (θ)e imϕ in Step 3 become coefficients of the form and the s fℓ −m (θ)e imϕ angular functions become coefficients of the form Since spin-weighted spherical harmonics are not orthogonal across different spins over the 2-sphere, one has to calculate these coefficients in general directly.Besides evaluating the integrals in Eqs. ( 130) and (131) for different (s 1 , ℓ 1 , m) and (s 2 , ℓ 2 , m) every time, there are also other approaches.One approach is to use the series-sum representation of s Y ℓm (θ, ϕ) in Eq. ( 61), as discussed in Appendix D with the results stored in a Mathematica notebook in [68].Now, the master equations of ϑ (1,1) and Ψ (1,1) 0,4 become completely radial, i.e., where s Hℓm is the radial Teukolsky operator for a spin s field in GR, and s A ℓm is the separation constant in the Teukolsky equation for a spin s field in GR, both of which can be found in [18].The coefficient s f k ℓm comes from the integral of s S ℓm (θ) and s f k ℓm (θ) over the 2-sphere, and similarly for its complex conjugate s fk ℓm .s f k ℓm and s fk ℓm will be given by Eqs. ( 130) and (131), respectively.
B. Radial part of the equation of ϑ (1,1)   In this subsection, we present the radial part of Eq. (36) in both the IRG and ORG found by following the procedures in Sec.IX A. In the IRG, we find the radial parts of terms proportional to e −iωt in Eqs. ( 93) and (95), respectively, to be where the terms proportional to b m ℓ,ℓ±1 or cos θ 0 Y ℓm (θ, ϕ) in Eq. ( 93) are at O(χ 2 ) after the angular integration.In the ORG, we find where s Rℓm (r) is the radial part of the Hertz potential given in Eq. ( 53), Λ ℓℓm 10s and Λ ℓℓm −10s are given by Eq. ( 130), and the prime denotes a derivative with respect to the radial coordinate r.The functions g ℓm i (r, ω, M ), h ℓm j (r, ω, M ), g ℓm i (r, ω, M ), h ℓm j (r, ω, M ) , where i ∈ [1,4] and j ∈ [1,2], are the same functions in Eqs. ( 93), ( 95), ( 112) and ( 113) and presented in a separate Mathematica notebook [68].
Using Eq. ( 57), we can replace the radial Hertz potential s Rℓm (r) and its derivative in Eqs. ( 133)-( 136) with s R (0,1) ℓm (r) and s R ′ (0, 1)  ℓm (r).Notice that the form of the equations remains similar with s Rℓm (r) now replaced by s R (0,1) ℓm (r) and the prefactors now new functions of {r, ω, M }.For instance, in the first parenthesis of Eq. ( 133), one finds that the prefactor of 2 R (0,1) ℓm (r) is (137) Each of the functions that would appear in V R ℓm (r), V □ ℓm (r), U R ℓm (r), and U □ ℓm (r) are separately presented in the supplementary Mathematica notebook due to their lengthy nature [68].
Combining Eq. ( 86) with Eqs. ( 133) and (134) [or Eqs. ( 135) and ( 136)], we now have a completely radial equation that describes the evolution of the scalar field perturbations, IRG: ORG: Recall that r s is the Schwarzschild radius, M is the mass of the BH, χ is the dimensionless spin parameter such that χ = a/M with a being the spin, 0 A ℓm is the separation constant for a spin-0 field [18], and V R ℓm (r), V □ ℓm (r), U R ℓm (r), U □ ℓm (r) are radial functions given in Eqs. ( 133 in the ORG] jointly.In [72], it was shown that one can turn Eq. ( 146) [or Eq. ( 150)] into an eigenvalue problem, following [24,38,39], such that the solutions of ηℓm correspond to the eigenvectors of the system, and the QNM frequencies ω ℓm are eigenvalues.
In the above equation, V †R ℓ −m refers to taking the complex conjugate of all the radial functions in s1s2c , Λ †ℓ1ℓ2m s1s2s }, and similarly for V †□ ℓ −m , U †R ℓ −m , and U †□ ℓ −m .Equations ( 138) and ( 139) can now be solved for to obtain the scalar-led QNM frequencies.Notice that there is a coupling between the scalar field perturbations and the gravitational perturbations in GR, which appear in the form of the Hertz potential radial function ±2 Rℓm (r) in the IRG or ORG, respectively.

C. Radial part of the equation of Ψ
(1,1) 0 In this subsection, we present the radial part of the modified Teukolsky equation for Ψ (1,1) 0 .Just like in the case of the ϑ field, the left-hand side of the modified Teukolsky equation Eq. ( 16) is the same as the Teukolsky equation in GR and is separable under the decomposition in Eq. (98), with its radial part given by Eqs. ( 96) and (97).To extract the radial part of the right-hand side, we follow the recipe provided in Sec.IX A to eliminate all angular dependence.
potential Ψ H in Eqs. ( 142) and ( 144) can be further expressed as functions of the radial Teukolsky function 2 R (0,1) ℓm (r) for the perturbed Ψ 0 in GR [18], as discussed in Sec.IX B. All necessary functions have been provided in a supplementary Mathematica notebook due to their lengthy nature [68].One can readily use existing wavefunction ansatz in the literature to evaluate the radial Teukolsky function [21,73].
Notice that the modified Teukolsky equation for the Weyl scalar perturbation Ψ (1,1) 0 exhibits coupling to the scalar field perturbation at O(ζ 1 , ϵ 1 ), but no such coupling is seen for the scalar field perturbation at the same order, unlike the case involving metric perturbations [9].
(147)-(149), we find As before, one has to solve the (ℓ, m) and (ℓ, −m) modes of Eq. ( 150) jointly to obtain the coefficient ηℓm in Eq. ( 128) and the QNM frequency ω ℓm .This analysis shows that the modified Teukolsky equations for the Weyl scalar perturbations Ψ (1,1) 0,4 and the scalar field perturbation ϑ (1,1) can be separated into a radial and angular piece.The radial piece can then be integrated numerically to obtain the QNM frequencies.

X. DISCUSSION
In this paper, we have employed the modified Teukolsky formalism in [23] to investigate the perturbations of slowly rotating BHs in dCS gravity at leading order in spin, where the BH spacetime is non-Ricci-flat, but remains of Petrov type D. To incorporate the slow-rotation approximation, we first extended the two-parameter expansion in [23] to a three-parameter expansion.Following [41,64], we then re-derived the null geodesics on the equatorial plane, from which we found the NP tetrad for slowly rotating BHs in dCS gravity up to O(χ).The resulting tetrad is the Kinnersly tetrad expanded to O(χ), with an additional adjustment accounting for the dCS correction.This tetrad is the same as the one in [64].Since BHs in dCS gravity are non-Ricci-flat, this direct extension of the Kinnersly tetrad leads to some nonzero background Weyl scalars Ψ 1 and Ψ 3 , so we performed additional tetrad rotations to remove them and computed all the background NP quantities in this rotated tetrad.
The source terms of the modified Teukolsky equation for Ψ 0,4 arise from two distinct contributions.Some of them originate from the homogeneous component of certain Bianchi and Ricci identities, so they only rely on the corrections to the background geometry.For Petrov type D spacetimes, these stationary corrections only couple to the perturbations of Ψ 0,4 , so we evaluated them using the NP quantities in the dCS background and the solutions to the Teukolsky equation in GR.The other source terms stem from the stress tensor associated with corrections to the Einstein-Hilbert action.In dCS gravity, these source terms couple the scalar field with the metric in GR.Thus, to completely evaluate them, we need to solve for the dynamical scalar field.In this case, we first evaluated the scalar field equation and used the same methodology to guide our calculations for Ψ 0,4 .
Since the scalar field is driven by dynamical metric perturbations in GR, one needs to first reconstruct the metric associated with curvature perturbations in GR.In this work, we chose to follow the CCK procedures developed in [33][34][35][36][59][60][61][62][63], where the perturbed metric is obtained from the Hertz potential, though other procedures in [37,41] may also apply.Using the reconstructed metric in [33][34][35][36][59][60][61][62][63], we then computed all the perturbed NP quantities in GR following the approach in [37,65].Since we also chose the gauge that the perturbations of Ψ 1,3 vanish in both GR and dCS gravity [23], we performed additional tetrad rotations to transform all the perturbed NP quantities into this gauge.In the end, projecting the scalar field equation onto the NP basis, we used the reconstructed NP quantities to express all the source terms as differential operators acting on the Hertz potential.
The Hertz potential can be obtained from the perturbations of Ψ 0,4 in GR, which are solutions to the Teukolsky equations in GR.Decomposing the Hertz potential into spin-weighted spheroidal harmonics, we presented the source terms of the scalar field equation in Boyer-Lindquist coordinates explicitly.The radial function of the Hertz potential was then determined from the radial function of the perturbed Ψ 0,4 in GR following [61].
To separate these master equations into radial and angular ordinary differential equations, we exploited the orthogonality properties of spin-weighted spheroidal harmonics and performed a harmonic decomposition to eliminate all angular dependence of the source terms.The homogeneous part of the scalar field equation naturally separates, so we obtain a purely radial differential equation [i.e., Eq. (138) in the IRG and Eq.(139) in the ORG].Similar procedures were then implemented for the modified Teukolsky equations of Ψ 0,4 .The source terms of the modified Teukolsky equations were expressed in terms of the Hertz potential and the dynamical scalar field.We then projected the source terms into the radial direction by integrating them over spin-weighted spheroidal harmonics.The homogeneous part of these equations separates in the same way as the Teukolsky equations in GR, so we also obtained two radial differential equations for Ψ 0 [i.e., Eq. ( 146)] and Ψ 4 [i.e., Eq. ( 150)], respectively.
Through these procedures, we obtained three coupled, ordinary (radial) differential equations for 2 R (1,1) ℓm (r) (or ℓm (r) in the ORG), where the first two are radials functions of Ψ (1,1) 0 ), and the last one is the radial function of ϑ (1,1) .All of these equations have the same structure.The left-hand side is the radial Teukolsky operator for particles of spin 2 (Ψ ), or spin 0 (ϑ (1,1) ).For the radial master equation of Ψ ), the right-hand side contains source terms that depend on Θ ℓ−m (r) ), where the last two are radial functions of Ψ (0,1) 0 (or ρ −4 Ψ (0,1) 4 ).For the radial master equation of ϑ (1,1) , the right-hand side contains source terms that depend on 2 R (0,1) ℓ−m (r) ).The coupled system therefore forms a (Sturm-Liouville) eigenvalue problem that should be amenable to standard procedures to find the eigenvectors and eigenvalues, i.e., the QNM and scalar frequencies.
The primary objective of this study was to apply the modified Teukolsky formalism in [23] to investigate perturbations of BHs in some specific modified theories of gravity.To illustrate this, we considered the case of slowly rotating BHs to leading order in spin within the framework of dCS gravity.Although the slow rotation approximation may not provide highly accurate results for more realistic BHs (with spins χ ∼ 0.6), it is a simplified problem for testing the newly developed formalism.Incorporating additional degrees of freedom associated with dCS gravity, coupled with the intricacies introduced by the metric reconstruction procedures, renders this calculation complex.Yet, in this work, we successfully demonstrated that the modified Teukolsky equation in [23] does not only decouple Weyl scalars Ψ 0,4 from other NP quantities but also admits a separation into radial and angular parts, a key advantage of the Teukolsky equation in GR, especially for rapidly rotating BHs.Although this paper focused on the first order in the slow rotation expansion, the separation of the modified Teukolsky equation should hold for any spin since the orthogonality properties of spin-weighted spheroidal harmonics we have used to separate the equation apply for a general spin.Thus, this calculation is an ideal initial step toward determining the QNM spectra for BHs with general spin in modified gravity.
This work creates a new path to directly calculate the corrections to the QNM frequencies for slowly rotating perturbed BHs in dCS gravity.Having obtained the master equations for the perturbed Weyl scalars Ψ 0,4 and the perturbed scalar field, we can now integrate these equations using numerical integration schemes, such as the eigenvalue perturbation method in [24,38,39] to find the QNM spectra.Moreover, the QNM spectra obtained using the modified Teukolsky formalism can then be compared to the results from the metric perturbation approach [6,7,9,31,74] and numerical relativity [75][76][77][78].Notice that, in higher-derivative gravity, Refs.[25,26] have followed our formalism to compute the QNMs in the slow-rotation expansion of BHs in that theory, and they obtained results valid for χ ≲ 0.7.Nonetheless, the dCS case we have focused on is more complicated due to the coupling to the scalar field equation.As discussed above, we have also presented in detail the angular dependence of the master equations of the perturbed Ψ 0,4 and ϑ and showed explicitly that they are separable, while Refs.[25,26] only briefly discussed using the orthogonality of spin-weighted spheroidal harmonics to extract the radial equations.
An additional aspect worth exploring is the phenomenon of isospectrality breaking in the QNM spectra.In GR, odd and even parity modes oscillate and decay at the same rate.However, certain modified theories of gravity have been shown to exhibit a breaking of isospectrality, e.g., dCS gravity [6,7,9,31,74], EdGB gravity [8,10,79], and higher-derivative gravity [55].The investigation of isospectrality breaking has, so far, primarily focused on metric perturbations, as the Zerilli-Moncrief and Regge-Wheeler equations naturally separate metric perturbations into even-and odd-parity sectors [14,15].However, for BHs with arbitrary spin, there are no known extensions of the Zerilli-Moncrief and the Regge-Wheeler equations, so we need to use the modified Teukolsky equation to study isospectrality breaking.In another study [72] involving all the authors, the definite-parity modes of curvature perturbations in modified gravity were found, and the features in these bGR theories that result in isospectrality breaking were revealed and demonstrated in several simple cases.Nonetheless, a direct mapping from the Zerilli-Moncrief and Regge-Wheeler functions to the modified Teukolsy equations of these definite-parity modes still remains unknown.The implementation of the modified Teukolsky equation in a concrete bGR theory has opened up possibilities for addressing these questions and more.
Building upon the insights gained from the present study, further investigations can be pursued involving more complex systems within various gravitational theories.As part of our collaborative effort, we are currently engaged in extending this calculation to derive the master equations and QNM spectra for BHs with arbitrary spin in dCS gravity, where the BH spacetime is Petrov type I.In addition, we are also actively involved in computing the master equations for rotating Petrov type I BHs within the framework of EdGB gravity.For the first time, we can explore the QNM spectra for BHs with general spin in a wide range of gravitational theories and spacetime geometries, which can then be compared with real observation data to scrutinize these possible deviations from GR.
For Weyl scalars at O(ζ 0 , ϵ 1 ), one can use Ricci identities to retrieve them from spin coefficients.The equations below work for both vacuum and non-vacuum spacetimes since we have linearly combined Ricci identities to remove NP Ricci scalars Φ ab following [37,65].
The equations above are precise, so one needs to linearize them when extracting the Weyl scalars at O(ζ 0 , ϵ 1 ) using the perturbed tetrad and spin coefficients at O(ζ 0 , ϵ 1 ).In Refs.[35,59], they also computed the perturbed Weyl scalars in the IRG and expressed them in terms of the Hertz potential, nn,1 = δ[1,3,0,−1] δ[0,4,0,3] ΨH .As we discussed in Sec.IV, our results using either the Ricci identities or direct computation from the linearized Riemann tensor agree for Ψ (0,1) 0,1,2,4 but not for Ψ (0,1) 3 due to different choices of the perturbed tetrad.Thus, we can use Eq.(B6) for Ψ (0,1) 0,1,2,4 and Eq.(B5d) for Ψ (0,1) 3 .For Schwarzschild, the equations above simplify into Ψ (0,1) 0 which reduce the number of independent components of C abcd to 10.These 10 independent components can be represented by 5 Weyl scalars and their conjugates.With complex conjugation, we can find all the components of the Weyl tensor using the symmetries above and the components below, where We have provided this series-sum representation of the coefficients in Eqs. ( 130) and (131) in a Mathematica notebook as Supplementary Material [68].

H
has the expansion in Eq. (

2
using our results in Sec.IV and Appendix B. The results of S

4 )
consists of the GR Teukolsky operator for a spin 2 (or spin −2) field acting on Ψ