Strong gravitational lensing occurs when photons pass through the vicinity of a black hole. We investigate this phenomenon in the context of a gravitationalwave event specifically when a black hole is settling into its final state. The deflection angle of photons mimics the ringdown pattern of the gravitational wave at intermediate times. At late times it has an inverse cubic dependence on observation time. The deviation angle increases exponentially as photons approach the photon ring orbit reflecting its unstable nature. Our findings are directly applicable to imaging scenarios involving stars against the background of compact binaries or circumbinary accretion disks particularly during the merger of two black holes.
['General Relativity and Quantum Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena']
Dynamical Lensing Tomography of Black Hole Ringdown
{'Dynamical Lensing Tomography of Black Hole Ringdown': "Zhen Zhong, 1, ∗ Vitor Cardoso, 2, 1, † and Yifan Chen 2, ‡  \n1  \nCENTRA, Departamento de F'ısica, Instituto Superior T'ecnico - IST, Universidade de Lisboa - UL, Avenida Rovisco Pais 1, 1049-001 Lisboa, Portugal 2 Niels Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  \nStrong gravitational lensing occurs when photons pass through the vicinity of a black hole. We investigate this phenomenon in the context of a gravitational-wave event, specifically when a black hole is settling into its final state. The deflection angle of photons mimics the ringdown pattern of the gravitational wave at intermediate times. At late times it has an inverse cubic dependence on observation time. The deviation angle increases exponentially as photons approach the photon ring orbit, reflecting its unstable nature. Our findings are directly applicable to imaging scenarios involving stars against the background of compact binaries, or circumbinary accretion disks, particularly during the merger of two black holes.", 'I. Introduction': "One of the foundational tests of General Relativity concerns the deflection of light by massive bodies [1-4]: the path of light rays is 'bent' by mass-energy, which distorts the spacetime fabric. In the last decade, it has become clear that the universe is flooded with traveling spacetime distortions, gravitational waves (GWs) traveling at the speed of light. Gravitational-wave astronomy provides a wealth of information on compact objects, particularly on the behavior of gravity in highly dynamical, high-curvature and redshift setups.  \nAn intriguing prospect concerns the detectability of light bending by GWs, which would add an extra information channel on strong gravity. This possibility was raised and studied by different authors in the weak field regime [5-9]. Our purpose here is to study light deflection in the strong-field, highly dynamical regime. In particular, we will study the appearance of a black hole (BH) which is approaching a quiescent state via relaxation in its characteristic modes of oscillation, called quasinormal modes (QNMs). The relaxation stage is also referred to as ringdown stage, and is a generic late-time description of the coalescence of compact objects [10, 11]. Our results could describe stars on the background of a binary BH merger, or the appearance of a circumbinary accretion disk surrounding two merging compact objects.  \nWe work in units where G = c = 1, and adopt the metric convention ( -, + , + , +).", 'II. Ray Tracing towards a Relaxing Black Hole': "A relaxing BH is well described by general relativistic first-order perturbation theory, where spacetime is expressed by small fluctuations away from the final, stationary state. Uniqueness results in General Relativity  \nsuggest that the final state, to a good approximation, is a Kerr BH [12-14]. We thus take our relaxing BH spacetime, of mass M and angular momentum aM to be given by  \ng µν = g (0) µν + h µν , (1)  \nwhere g (0) µν is the Kerr metric (we will use Boyer-Lindquist ( t, r, θ, ϕ ) coordinates throughout [15]). The term h µν denotes first-order perturbations.  \nInstead of working directly with the metric perturbations h µν , we consider instead curvature perturbations [16, 17]. These can be separated and decoupled and all the relevant information about radiative degrees of freedom is encapsulated in a master variable  \nψ s ( t, r, θ, ϕ ) = s Z √ 2 π s S ( θ ) s R ( r ) e -iωt + imϕ . (2)  \nThe gravitational quantity of interest is Ψ (1) 4 = ψ -2 ( t, r, θ, ϕ )( r -ia cos θ ) -4 , the first-order perturbation of one of the Weyl scalars, which is directly related to the amplitude of GWs far from the source. Notation and conventions are shown in Supplemental Material A 1. The function s S ( θ ) is a spin-weighted spheroidal harmonic, while the radial Teukolsky function s R ( r ) obeys a homogeneous second-order differential equation. We fix its normalization by insisting on the large-distance behavior. Then, s Z controls the amplitude of the GW, by matching the ringdown stage to some physically generated process. The characteristic modes are found by imposing proper boundary conditions on s R ( r ) , which selects only a discrete set of frequencies ω , the ringdown frequencies [10]. All the functions and quantities carry a dependence on angular numbers ℓ, m and a possible overtone number n [10, 16, 17]. For simplicity, we do not include these indices in the subscripts. Although a true GW signal comprises a superposition of various modes, our discussion focuses on the dominant ℓ = m = 2 mode. Additionally, we exclusively address the n = 0 mode, which lasts longer compared to the more rapidly decaying higher overtones. For QNMs characterized by frequency ω and labeled by ( ℓ, m, n ), there exists a corresponding  \n'mirror' mode with frequency -¯ ω and ( ℓ, -m,n ). In the final spin frame, these paired modes primarily represent GW emissions propagating in opposite directions along the spin axis-one mode predominantly radiating towards the northern hemisphere and the other towards the southern hemisphere [18-21]. The observed waveform typically manifests as a superposition of these two distinct sets of modes. However, one of the paired modes may be less prominent due to weaker excitation [19, 22]. In the following analysis, we primarily focus on the m = 2 mode with Re( ω ) > 0, which predominantly emits in the southern hemisphere, as more general cases exhibit similar features.  \nOnce -2 R ( r ) is known, we reconstruct the metric fluctuation h µν following well-known prescriptions [23-25], summarized in Supplemental Material A 2. The QNMs diverge as h µν ∼ e Im( ω )( t -r ∗ ) when t -r ∗ →-∞ , where Im( ω ) represents the imaginary part of the ringdown frequency ω and r ∗ is the tortoise coordinate [10]. To mitigate the divergence, we truncate the perturbation metric h µν by applying a Heaviside-like filter function H ( t -r ∗ ):  \nH ( t -r ∗ ) = 1 2 ( tanh t -r ∗ q +1 ) , (3)  \nto obtain H ( t -r ∗ ) h µν . Here, 0 < q < -2 / Im( ω ), where q controls the intermediate region width of the filter. Although various values of q were tested, the modifications had only a minor impact on the outcome, especially regarding the late-time frequency behavior. In the following, we employ a relatively sharp truncation in the filter function H ( t -r ∗ ), with -q Im( ω ) = 0 . 1.  \nDuring the ringdown phase, a photon trajectory is influenced by the metric perturbations encountered along its path, governed by the geodesic equation:  \nd 2 x µ d τ 2 +Γ µ νρ d x ν d τ d x ρ d τ = 0 , (4)  \nwhere Γ denotes the Levi-Civita connection, and τ represents the proper time.  \nTo characterize the observable effects, it is essential to define the observer frame. In an unperturbed Kerr BH scenario, a zero angular momentum observer (ZAMO) frame is typically used [26]. However, in a dynamic spacetime, the ZAMO tetrads are neither orthogonal nor normalized. To rectify this, we apply a Gram-Schmidt orthogonalization process to establish a normalized, orthogonal tetrad basis derived from the ZAMO tetrads, as detailed in Supplemental Material B 1. Using this newly defined observer frame, we initiate geodesics towards the BH employing backward ray tracing. Here, each pixel on the observer plane, represented by coordinates ( ˜ X, ˜ Y ) [27-29], corresponds to photon momentum in this frame. The photon's 4-momentum in the dynamic spacetime p µ ≡ d x µ / d τ is then constructed from this orthogonal tetrad basis, similar to the procedure established in Kerr spacetime [30]. For convenience, we also introduce polar coordinates for the observer plane  \ncentered on the direction directly facing the BH in Kerr spacetime, denoted as (˜ ρ, ˜ φ ), where ˜ ρ is the impact parameter and ˜ φ is the angle.  \nPhoton geodesics in the relaxing Kerr spacetime are then computed. We position an observer at r o = 10 9 M , confirming convergence of results at large r o due to the decay of h µν amplitude as 1 /r at greater distances. Notice that even though we employ backward ray tracing, photon propagation in the + t direction still aims towards the observer, aligning with the GWs that co-propagate in the same direction upon departure from the BH.", 'III. Lensing Tomography': "We explore geodesics in a relaxing Kerr spacetime, particularly focusing on the merger events of supermassive binary BHs (SMBHBs) with equal mass ratios. The remnant BH from each mass ratio merger events typically characterized by an angular momentum parameter a ∼ 0 . 7 M [31], a value consistent with numerous observations [32, 33] and considered as our fiducial value. The GW signal from these mergers is dominated by the ℓ = | m | = 2 mode [10, 19, 34, 35], with the n = 0 mode persisting the longest. This fundamental mode is characterized by Mω = 0 . 5326 -0 . 08079 i [10, 36].  \nIn the vicinity of the BH, the metric perturbation amplitude can reach up to M 2 Ψ (1) 4 ∼ 0 . 1 [34]. The total radiated energy during the ringdown is limited to approximately 3% of the remnant BH's mass [34], setting an upper limit on the magnitude of | -2 Z | for the ℓ = m = 2 mode at about 0 . 24. We adopt this upper bound as the fiducial value for our analyses, our results can be trivially re-scaled for generic amplitudes.  \nWe first consider a face-on observer with θ o ≈ 0. The top panel of Fig. 1 illustrates backward photon trajectories, all launched from the same direction, using identical backward impact parameters relative to the observer positioned at r 0 = 10 9 M . These trajectories are depicted using the Cartesian Kerr-Schild (KS) coordinate system ( x, y, z ). The black line represents geodesics in an unperturbed Kerr spacetime, while various colors indicate different emission times t o -r ∗ ( r o ). The bottom panel displays the growth of deviations between unperturbed and deflected geodesics, defined as ∆ KS ≡ (∆ x 2 +∆ y 2 +∆ z 2 ) 1 / 2 , representing the distance between the two trajectories at the same backward propagating time t 0 -t -r ∗ ( r o ), colors corresponding to trajectory lines in the top panel.  \nIn Fig. 1, we select ˜ ρ = 5 . 041 M as the (backward) impact parameter, which is close to the critical curve ˜ ρ c ≈ 5 . 04 M for this BH spin. This curve, defined in Kerr spacetime, represents the geodesics that can propagate around the BH indefinitely. Some of the perturbed geodesics, shown in dashed lines, terminate on the BH instead of propagating towards infinity, indicating a distortion of the critical curve during the ringdown.  \nThe evolution of geodesic deviations, launched at dif-  \nFigure 2 displays the deflected angles in polar (∆ θ ) and azimuthal (∆ ϕ ) components as a function of observer frame time, t 0 -r ∗ ( r o ), for trajectories initially launched from various polar coordinate angles ˜ φ on the observer plane. The impact parameters are fixed at ˜ ρ = 5 . 4 M . Due to ringdown metric perturbations, axial symmetry is not preserved even in face-on observations. Notably, both ∆ ϕ and ∆ θ show identical variations for values of ˜ φ separated by π . This pattern is precisely attributable to the m = 2 mode of the metric perturbations, akin to those  \n<!-- image -->  \n10  \n10  \n-  \n10  \n-  \n10  \n-  \n1  \n1  \n3  \n5  \nM  \n/  \nKS  \n∆  \n0  \n50  \n100  \n-  \n(  \nt  \no  \n-  \nt  \nr  \n∗  \n(  \nr  \no  \n))  \n/M  \nFIG. 1: Top: Illustration of photon trajectories in a relaxing Kerr BH spacetime with spin a = 0 . 7 M , viewed from a face-on observer with θ 0 ≈ 0. The metric perturbation is modeled using the n = 0 and ℓ = m = 2 mode, with magnitude normalized to | -2 Z | = 0 . 24. The initial impact parameters are set at ˜ ρ = 5 . 041 M , slightly exceeding the critical curve ˜ ρ c = 5 . 04 M in Kerr spacetime. The black line represents unperturbed geodesics in Kerr spacetime, an arrow indicating the direction from the observer to the BH to depict the backward ray tracing method used. Different colors represent geodesics launched at various observation times t 0 -r ∗ ( r 0 ), detailed in the bottom panel. The dashed line marks a trajectory terminating on the BH, whereas the others extend towards infinity. The gray plane indicates the BH equatorial plane. Bottom: Growth of geodesic deviation ∆ KS representing the distance between unperturbed and deflected geodesics at the same backward propagating time t 0 -t -r ∗ ( r 0 ), shown in Cartesian KS coordinates.  \n(  \nt  \no  \n-  \nr  \n∗  \n(  \nr  \no  \n))  \n/M  \n0  \n10  \n20  \n30  \n40  \n50  \n<!-- image -->  \n-  \n∗  \nFIG. 2: Variations of the deflected angles in the azimuthal (∆ ϕ ) and polar (∆ θ ) components for geodesics terminating at infinity, plotted as a function of observation time t 0 -r ∗ ( r o ). The BH, observer, and ringdown configuration follow the setup described in Fig. 1. Different colors represent geodesics launched from various angles ˜ φ on the polar coordinate of the observer plane, with impact parameters fixed at ˜ ρ = 5 . 4 M . Variations for cases separated by π in ˜ φ appear identical, reflecting the m = 2 nature of the BH ringdown.  \nferent observer times, initially shows a weak increase after departing from the observer. Geodesics launched later exhibit deviations exponentially smaller than those launched earlier, reflecting the exponential decay of GW amplitude near the BH during the ringdown phase. Closer to the BH, all deviations begin to exhibit exponential growth due to the instability of the radial potential in the light ring region [29, 37-47]. Here, the interplay between light ring instability and GW waves induces addiontial oscillation features in the exponential growth. The geodesics that depart from the light ring region begin to develop deviations that increase linearly with time.  \nTo quantify the lensing during ringdown, we analyze backward photons that ultimately reach infinity, corresponding to far-away sources. We assess the angular difference between the asymptotic direction for Kerr geodesic and the deflected geodesic as observables, including the grand arc difference ˜ δ , and variations in the azimuthal angle ∆ ϕ and polar angle ∆ θ .  \ngenerated by vector or tensor superradiant clouds [47], demonstrating that the deflected angles can effectively dissect the polarization pattern of the GWs around the BH.  \nNext, we adjust the impact parameters ∆˜ ρ ≡ ˜ ρ -˜ ρ c relative to the critical curve ˜ ρ c ≈ 5 . 04 M , as defined in Kerr spacetime. In Fig. 3, we illustrate variations of the deflected angle ˜ δ . Due to the implementation of the filter, deflection begins around t 0 -r ∗ ( r o ) ≈ 0. For impact parameters close to the critical curve, specifically in cases where ∆˜ ρ ≤ 1 . 5 M , the early time evolution of the deflected angle exhibits a pattern of exponential decay and oscillation. The timescales of this decay and the oscillation period align with the inverses of the imaginary and real parts of ω , respectively, suggesting that these observables can effectively capture key characteristics of the BH ringdown.  \nFurthermore, when comparing the overall amplitude of ˜ δ , there is a noticeable exponential increase dependent on ∆˜ ρ . This enhancement is again attributed to the instability of the BH light ring region, which causes an exponential growth in the perturbed deviation during propagation around the BH [29, 37-45, 47]. As ˜ ρ approaches the critical curve, this instability intensifies. The instability is quantified by the square root of the second derivative of the radial potential [45], which increases linearly as ∆˜ ρ decreases, thereby leading to an exponential rise in the deflected angle.  \nOn the other hand, for late behavior and cases with large impact parameters ∆˜ ρ/M = 5 , 10 , 15, the variations exhibit a slow decay without oscillations. The two different types of behavior is attributable to the deflections developing in two distinct regions. From the BH to the observer in the + t direction, the GW and photons co-propagate, so the deflection that develops along this path reflects the ringdown metric perturbations present when the photon starts to co-propagate with the GWs. Before entering the co-propagation region, particularly as the photon travels from infinity towards the BH, the deflection received along the path is predominantly influenced by the time it passes through the filter at t ≈ r ∗ . For smaller impact parameters, the co-propagation duration lasts longer, resulting in deflected angles starting with an exponential decay and oscillation, mirroring the GW ringdown. However, in late-time scenarios and with larger impact parameters, the deflection is primarily influenced by non-co-propagating impacts after t ≈ r ∗ . The case with ∆˜ ρ = 15 M shows larger deflection than those with ∆˜ ρ = 5 M and 10 M because its polar angle θ at t ≈ r ∗ is closer to π , where the metric perturbation for the ℓ = m = 2 mode is most dominant. Our results are consistent with a power law, ˜ δ ∼ ( t o -r ∗ ) -p , where p ≈ 3. This scaling behavior resembles that observed in GW lensing from localized sources [9, 48-50], where the dominant contribution to lensing at distances far from the source is proportional to h µν p µ p ν , which is only non-vanishing for the O (1 /r 3 ) component of the metric perturbations.  \nFIG. 4: Similar to Fig. 3, this figure displays variations in the deflected angles of geodesics terminating at infinity, as observed from different inclination angles θ 0 and plotted against observation time t 0 -r ∗ ( r o ). The impact parameter is fixed at ˜ ρ = 10 M on the + ˜ X axis.  \n<!-- image -->  \n-  \n∗  \nFIG. 3: Deflection angle of geodesics terminating at infinity, defined as the grand arc difference with the unperturbed direction, plotted as a function of the observation time t 0 -r ∗ ( r o ). Different colors represent geodesics launched with various impact parameter shifts ∆˜ ρ ≡ ˜ ρ -˜ ρ c from the critical curve in Kerr spacetime, where ˜ ρ c ≈ 5 . 04 M , along the + ˜ X -axis. The BH, observer, and ringdown configuration are consistent with the setup described in Fig. 1. Near the critical curve, we find an oscillatory, exponentially decaying dependence at intermediate times, in very good (to sub-percent level) agreement with the ringdown of the GW itself. At later times, a pattern of slow decay is observed, consistent with a power law ˜ δ ∼ ( t o -r ∗ ) -p with p ∼ 3.  \n<!-- image -->  \n-  \n∗  \nWe also display cases with different observer inclination angles θ o in Fig. 4, fixing the impact parameters at  \n˜ ρ = 10 M along the + ˜ X -axis. Similar to the observations in Fig. 3, the variations exhibit two types of behavior: a slow decay predominates at smaller θ o , while larger θ o are primarily characterized by ringdown-like decay and oscillation. This distinction can again be attributed to the propagation of the polar angle from infinity to the BH, where smaller θ o corresponds to larger θ in this trajectory. For θ o = π (1 -10 -4 ) for example, the deflection angle mimics the ringdown frequency to one part in 10 10 !  \nFurthermore, we consider the effect of incorporating the QNM's mirror mode with equal energy, reflecting a more realistic scenario where the SMBHB was in a quasicircular orbit before the merger. This setup also exhibits the two previously discussed characteristics: an initial rapid decay during the ringdown phase, followed by a late-time power-law tail. Notably, the ringdown phase becomes more pronounced as the metric perturbations are now comparable at both θ ∼ 0 and θ ∼ π .", 'IV. Discussion': "This study has explored the dynamic lensing of photons during the ringdown phase following the merger of SMBHBs. By calculating photon trajectories around a relaxing Kerr BH using backward ray tracing, we demonstrate how metric perturbations during the ringdown significantly influence photon paths. In some instances, variations in photon deflection angles closely correspond to the GW ringdown behavior, particularly when photons co-propagate with GWs from the BH to the observer, thereby offering insights into the fascinating dynamics of spacetime.  \nThe visualization of ringdown through photon lensing opens a new avenue for multi-messenger astronomy. While pulsar timing arrays [51-54], astrometry [55-57], binary orbits [58, 59], and fast radio burst timing [60] may directly detect ringdown GWs from SMBHB mergers, observations from the Event Horizon Telescope [61, 62] and future upgrades [63], including space missions [64-66], could provide visualizations of a 'Waltz' between GW and photon co-propagation. Each pixel on the observer plane could reveal additional insights by dissecting the metric perturbations in the vicin-  \n- [1] Clifford M. Will, 'The Confrontation between General Relativity and Experiment,' Living Rev. Rel. 17 , 4 (2014), arXiv:1403.7377 [gr-qc].\n- [2] Jos'e P. S. Lemos, Carlos A. R. Herdeiro, and Vitor Cardoso, 'Einstein and Eddington and the eclipse in Principe: Celebration and science 100 years after,' (2019), arXiv:1912.08354 [physics.hist-ph].\n- [3] Matthias Bartelmann and Peter Schneider, 'Weak gravitational lensing,' Phys. Rept. 340 , 291-472 (2001), arXiv:astro-ph/9912508.\n- [4] Matthias Bartelmann, 'Gravitational Lensing,' Class.  \nity of the relaxing BH.  \nClearly, even though we focus on ringdown, their quasiperiodic nature indicates that our results can be taken as indicative of what happens during earlier stages in the life of a binary emitting gravitational radiation.  \nWhile our use of backward ray tracing sometimes leads to the disappearance or replacement of a new source at a given pixel, predicting the apparent position shift of a source remains challenging. A future research direction will involve studying photon propagation from a given source point towards the observer in a dynamic spacetime. This endeavor will require a novel approach termed forward ray tracing, presently only feasible by leveraging conserved quantities in Kerr spacetime [67].", 'Acknowledgments': "Wethank Maarten van de Meent for helpful discussions and for providing us with useful routines for metric reconstruction. We also extend our thanks to Michael Boyle for his insightful discussions on spin-weighted spherical harmonics. Z.Z. acknowledges financial support from China Scholarship Council (No. 202106040037). V.C. is a Villum Investigator and a DNRF Chair. We acknowledge support by VILLUM Foundation (grant no. VIL37766) and the DNRF Chair program (grant no. DNRF162) by the Danish National Research Foundation. We acknowledge financial support provided under the European Union's H2020 ERC Advanced Grant 'Black holes: gravitational engines of discovery' grant agreement no. Gravitas-101052587. Views and opinions expressed are however those of the author only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 101007855 and No 101131233. Y.C. acknowledges the support of the Rosenfeld foundation and the European Consortium for Astroparticle Theory in the form of an Exchange Travel Grant.  \nQuant. Grav. 27 , 233001 (2010), arXiv:1010.3829 [astroph.CO].  \n- [5] A. Labeyrie, 'Lensing effects of gravitational radiation near celestial sources.' Astronomy and Astrophysics 268 , 823-828 (1993).\n- [6] Redouane Fakir, 'Gravity Wave Watching,' Astrophys. J. 426 , 74 (1994), arXiv:gr-qc/9304003 [gr-qc].\n- [7] Ruth Durrer, 'Light deflection in perturbed friedmann universes,' Phys. Rev. Lett. 72 , 3301-3304 (1994).\n- [8] Nick Kaiser and Andrew H. Jaffe, 'Bending of light by gravity waves,' Astrophys. J. 484 , 545-554 (1997),\n- arXiv:astro-ph/9609043.\n- [9] Thibault Damour and Gilles Esposito-Farese, 'Light deflection by gravitational waves from localized sources,' Phys. Rev. D 58 , 044003 (1998), arXiv:gr-qc/9802019.\n- [10] Emanuele Berti, Vitor Cardoso, and Andrei O. Starinets, 'Quasinormal modes of black holes and black branes,' Class. Quant. Grav. 26 , 163001 (2009), arXiv:0905.2975 [gr-qc].\n- [11] Vishal Baibhav, Mark Ho-Yeuk Cheung, Emanuele Berti, Vitor Cardoso, Gregorio Carullo, Roberto Cotesta, Walter Del Pozzo, and Francisco Duque, 'Agnostic black hole spectroscopy: Quasinormal mode content of numerical relativity waveforms and limits of validity of linear perturbation theory,' Phys. Rev. D 108 , 104020 (2023), arXiv:2302.03050 [gr-qc].\n- [12] Piotr T. Chrusciel, Joao Lopes Costa, and Markus Heusler, 'Stationary Black Holes: Uniqueness and Beyond,' Living Rev. Rel. 15 , 7 (2012), arXiv:1205.6112 [gr-qc].\n- [13] Vitor Cardoso and Leonardo Gualtieri, 'Testing the black hole 'no-hair' hypothesis,' Class. Quant. Grav. 33 , 174001 (2016), arXiv:1607.03133 [gr-qc].\n- [14] Vitor Cardoso, Kyriakos Destounis, Francisco Duque, Rodrigo Panosso Macedo, and Andrea Maselli, 'Black holes in galaxies: Environmental impact on gravitationalwave generation and propagation,' Phys. Rev. D 105 , L061501 (2022), arXiv:2109.00005 [gr-qc].\n- [15] James M. Bardeen, William H. Press, and Saul A Teukolsky, 'Rotating black holes: Locally nonrotating frames, energy extraction, and scalar synchrotron radiation,' Astrophys. J. 178 , 347 (1972).\n- [16] Saul A. Teukolsky, 'Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electromagnetic and neutrino field perturbations,' Astrophys. J. 185 , 635-647 (1973).\n- [17] S. A. Teukolsky, 'Rotating black holes - separable wave equations for gravitational and electromagnetic perturbations,' Phys. Rev. Lett. 29 , 1114-1118 (1972).\n- [18] Xiang Li, Ling Sun, Rico Ka Lok Lo, Ethan Payne, and Yanbei Chen, 'Angular emission patterns of remnant black holes,' Phys. Rev. D 105 , 024016 (2022), arXiv:2110.03116 [gr-qc].\n- [19] Emanuele Berti, Vitor Cardoso, and Clifford M. Will, 'On gravitational-wave spectroscopy of massive black holes with the space interferometer LISA,' Phys. Rev. D 73 , 064030 (2006), arXiv:gr-qc/0512160.\n- [20] Gregory B. Cook and Maxim Zalutskiy, 'Gravitational perturbations of the Kerr geometry: High-accuracy study,' Phys. Rev. D 90 , 124021 (2014), arXiv:1410.7698 [gr-qc].\n- [21] William Krivan, Pablo Laguna, Philippos Papadopoulos, and Nils Andersson, 'Dynamics of perturbations of rotating black holes,' Phys. Rev. D 56 , 3395-3404 (1997), arXiv:gr-qc/9702048.\n- [22] Eleanor Hamilton, Lionel London, and Mark Hannam, 'Ringdown frequencies in black holes formed from precessing black-hole binaries,' Phys. Rev. D 107 , 104035 (2023), arXiv:2301.06558 [gr-qc].\n- [23] P. L. Chrzanowski, 'Vector Potential and Metric Perturbations of a Rotating Black Hole,' Phys. Rev. D 11 , 2042-2062 (1975).\n- [24] J. M. Cohen and L. S. Kegeles, 'Electromagnetic fields in curved spaces - a constructive procedure,' Phys. Rev. D 10 , 1070-1084 (1974).  \n- [25] L. S. Kegeles and J. M. Cohen, 'CONSTRUCTIVE PROCEDURE FOR PERTURBATIONS OF SPACETIMES,' Phys. Rev. D 19 , 1641-1664 (1979).\n- [26] V. P. Frolov and I. D. Novikov, eds., Black hole physics: Basic concepts and new developments (1998).\n- [27] J. M. Bardeen, 'Timelike and null geodesics in the Kerr metric,' Proceedings, Ecole d'Et'e de Physique Th'eorique: Les Astres Occlus : Les Houches, France, August, 1972, 215-240 , 215-240 (1973).\n- [28] Samuel E. Gralla, Alexandru Lupsasca, and Andrew Strominger, 'Observational Signature of High Spin at the Event Horizon Telescope,' Mon. Not. Roy. Astron. Soc. 475 , 3829-3853 (2018), arXiv:1710.11112 [astro-ph.HE].\n- [29] Samuel E. Gralla and Alexandru Lupsasca, 'Lensing by Kerr Black Holes,' Phys. Rev. D 101 , 044031 (2020), arXiv:1910.12873 [gr-qc].\n- [30] Pedro V. P. Cunha, Carlos A. R. Herdeiro, Eugen Radu, and Helgi F. Runarsson, 'Shadows of Kerr black holes with and without scalar hair,' Int. J. Mod. Phys. D 25 , 1641021 (2016), arXiv:1605.08293 [gr-qc].\n- [31] Alessandra Buonanno, Gregory B. Cook, and Frans Pretorius, 'Inspiral, merger and ring-down of equal-mass black-hole binaries,' Phys. Rev. D 75 , 124018 (2007), arXiv:gr-qc/0610122.\n- [32] Alexander H. Nitz, Collin Capano, Alex B. Nielsen, Steven Reyes, Rebecca White, Duncan A. Brown, and Badri Krishnan, '1-OGC: The first open gravitationalwave catalog of binary mergers from analysis of public Advanced LIGO data,' Astrophys. J. 872 , 195 (2019), arXiv:1811.01921 [gr-qc].\n- [33] B. P. Abbott et al. (LIGO Scientific, Virgo), 'GWTC1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs,' Phys. Rev. X 9 , 031040 (2019), arXiv:1811.12907 [astro-ph.HE].\n- [34] Emanuele Berti, Vitor Cardoso, Jose A. Gonzalez, Ulrich Sperhake, Mark Hannam, Sascha Husa, and Bernd Bruegmann, 'Inspiral, merger and ringdown of unequal mass black hole binaries: A Multipolar analysis,' Phys. Rev. D 76 , 064034 (2007), arXiv:gr-qc/0703053.\n- [35] Michele Maggiore, Gravitational Waves: Volume 2: Astrophysics and Cosmology (Oxford University Press, 2018).\n- [36] Leo C. Stein, 'qnm: A Python package for calculating Kerr quasinormal modes, separation constants, and spherical-spheroidal mixing coefficients,' J. Open Source Softw. 4 , 1683 (2019), arXiv:1908.10377 [gr-qc].\n- [37] M. A. Podurets, 'Asymptotic Behavior of the Optical Luminosity of a Star in Gravitational Collapse,' Soviet Ast. 8 , 868 (1965).\n- [38] William L. Ames and Kip S. Thorne, 'The Optical Appearance of a Star that is Collapsing Through its Gravitational Radius,' Astrophys. J. 151 , 659 (1968).\n- [39] J.-P. Luminet, 'Image of a spherical black hole with thin accretion disk,' Astron. Astrophys. 75 , 228-235 (1979).\n- [40] G.A. Campbell and R.A. Matzner, 'A model for peaking of galactic gravitational radiation,' J. Math. Phys. 14 , 1-6 (1973).\n- [41] C. T. Cunningham and J. M. Bardeen, 'The Optical Appearance of a Star Orbiting an Extreme Kerr Black Hole,' Astrophys. J. Lett. 173 , L137 (1972).\n- [42] Heino Falcke, Fulvio Melia, and Eric Agol, 'Viewing the shadow of the black hole at the galactic center,' Astrophys. J. Lett. 528 , L13 (2000), arXiv:astro-ph/9912263.\n- [43] S. Chandrasekhar, The Mathematical Theory of Black Holes (Oxford University Press, New York, 1983).\n- [44] Vitor Cardoso, Alex S. Miranda, Emanuele Berti, Helvi Witek, and Vilson T. Zanchin, 'Geodesic stability, Lyapunov exponents and quasinormal modes,' Phys. Rev. D 79 , 064016 (2009), arXiv:0812.1806 [hep-th].\n- [45] Michael D. Johnson et al. , 'Universal interferometric signatures of a black hole's photon ring,' Sci. Adv. 6 , eaaz1310 (2020), arXiv:1907.04329 [astro-ph.IM].\n- [46] Vitor Cardoso, Francisco Duque, and Arianna Foschi, 'Light ring and the appearance of matter accreted by black holes,' Phys. Rev. D 103 , 104044 (2021), arXiv:2102.07784 [gr-qc].\n- [47] Yifan Chen, Xiao Xue, Richard Brito, and Vitor Cardoso, 'Photon Ring Astrometry for Superradiant Clouds,' Phys. Rev. Lett. 130 , 111401 (2023), arXiv:2211.03794 [gr-qc].\n- [48] Sergei M. Kopeikin, Gerhard Schaefer, Carl R. Gwinn, and T. Marshall Eubanks, 'Astrometric and timing effects of gravitational waves from localized sources,' Phys. Rev. D 59 , 084023 (1999), arXiv:gr-qc/9811003.\n- [49] Sergei M. Kopeikin and Gerhard Schaefer, 'Lorentz covariant theory of light propagation in gravitational fields of arbitrary moving bodies,' Phys. Rev. D 60 , 124002 (1999), arXiv:gr-qc/9902030.\n- [50] Sergei Kopeikin, Pavel Korobkov, and Alexander Polnarev, 'Propagation of light in the field of stationary and radiative gravitational multipoles,' Class. Quant. Grav. 23 , 4299-4322 (2006), arXiv:gr-qc/0603064.\n- [51] Gabriella Agazie et al. (NANOGrav), 'The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background,' Astrophys. J. Lett. 951 , L8 (2023), arXiv:2306.16213 [astro-ph.HE].\n- [52] J. Antoniadis et al. (EPTA, InPTA:), 'The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals,' Astron. Astrophys. 678 , A50 (2023), arXiv:2306.16214 [astro-ph.HE].\n- [53] Daniel J. Reardon et al. , 'Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,' Astrophys. J. Lett. 951 , L6 (2023), arXiv:2306.16215 [astro-ph.HE].\n- [54] Heng Xu et al. , 'Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I,' Res. Astron. Astrophys. 23 , 075024 (2023), arXiv:2306.16216 [astro-ph.HE].\n- [55] V. B. Braginsky, N. S. Kardashev, I. D. Novikov, and A. G. Polnarev, 'Propagation of electromagnetic radiation in a random field of gravitational waves and space radio interferometry,' Nuovo Cim. B 105 , 1141-1158 (1990).\n- [56] Laura G. Book and Eanna E. Flanagan, 'Astrometric Effects of a Stochastic Gravitational Wave Background,' Phys. Rev. D 83 , 024024 (2011), arXiv:1009.4192 [astroph.CO].\n- [57] Christopher J. Moore, Deyan P. Mihaylov, Anthony Lasenby, and Gerard Gilmore, 'Astrometric Search Method for Individually Resolvable Gravitational Wave Sources with Gaia,' Phys. Rev. Lett. 119 , 261102 (2017), arXiv:1707.06239 [astro-ph.IM].\n- [58] Diego Blas and Alexander C. Jenkins, 'Detecting stochastic gravitational waves with binary resonance,' Phys. Rev. D 105 , 064021 (2022), arXiv:2107.04063 [grqc].\n- [59] Diego Blas and Alexander C. Jenkins, 'Bridging the  \n- µ Hz Gap in the Gravitational-Wave Landscape with Binary Resonances,' Phys. Rev. Lett. 128 , 101103 (2022), arXiv:2107.04601 [astro-ph.CO].\n- [60] Zhiyao Lu, Lian-Tao Wang, and Huangyu Xiao, 'A New Probe of µ Hz Gravitational Waves with FRB Timing,' (2024), arXiv:2407.12920 [gr-qc].\n- [61] Kazunori Akiyama et al. (Event Horizon Telescope), 'First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole,' Astrophys. J. Lett. 875 , L1 (2019), arXiv:1906.11238 [astro-ph.GA].\n- [62] Kazunori Akiyama et al. (Event Horizon Telescope), 'First Sagittarius A* Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole in the Center of the Milky Way,' Astrophys. J. Lett. 930 , L12 (2022), arXiv:2311.08680 [astro-ph.HE].\n- [63] D. Ayzenberg et al. , 'Fundamental Physics Opportunities with the Next-Generation Event Horizon Telescope,' (2023), arXiv:2312.02130 [astro-ph.HE].\n- [64] Kari Haworth et al. , 'Studying black holes on horizon scales with space-VLBI,' (2019), arXiv:1909.01405 [astro-ph.IM].\n- [65] Michael D. Johnson et al. , 'The Black Hole Explorer: Motivation and Vision,' (2024) arXiv:2406.12917 [astroph.IM].\n- [66] Kazunori Akiyama et al. , 'The Japanese Vision for the Black Hole Explorer Mission,' (2024) arXiv:2406.09516 [astro-ph.IM].\n- [67] Lihang Zhou, Zhen Zhong, Yifan Chen, and Vitor Cardoso, 'Forward Ray Tracing and Hot Spots in Kerr Spacetime,' (2024), arXiv:2408.16049 [gr-qc].\n- [68] Maarten van de Meent, 'Gravitational self-force on eccentric equatorial orbits around a Kerr black hole,' Phys. Rev. D 94 , 044034 (2016), arXiv:1606.06297 [gr-qc].\n- [69] Andrew Spiers, Adam Pound, and Jordan Moxon, 'Second-order Teukolsky formalism in Kerr spacetime: formulation and nonlinear source,' (2023), arXiv:2305.19332 [gr-qc].\n- [70] Adam Pound and Barry Wardell, 'Black hole perturbation theory and gravitational self-force,' in Handbook of Gravitational Wave Astronomy , edited by Cosimo Bambi, Stavros Katsanevas, and Konstantinos D. Kokkotas (Springer Singapore, Singapore, 2020) pp. 1119.\n- [71] Robert M. Wald, 'Construction of Solutions of Gravitational, Electromagnetic, Or Other Perturbation Equations from Solutions of Decoupled Equations,' Phys. Rev. Lett. 41 , 203-206 (1978).\n- [72] Oscar J. C. Dias, Harvey S. Reall, and Jorge E. Santos, 'Kerr-CFT and gravitational perturbations,' JHEP 08 , 101 (2009), arXiv:0906.2380 [hep-th].\n- [73] R. Cartas-Fuentevilla, 'Perturbations of EinsteinMaxwell dilaton fields,' Phys. Rev. D 56 , 7700-7707 (1997).\n- [74] Vahid Toomani, Peter Zimmerman, Andrew Spiers, Stefan Hollands, Adam Pound, and Stephen R. Green, 'New metric reconstruction scheme for gravitational self-force calculations,' Class. Quant. Grav. 39 , 015019 (2022), arXiv:2108.04273 [gr-qc].\n- [75] Stephen R. Green, Stefan Hollands, and Peter Zimmerman, 'Teukolsky formalism for nonlinear Kerr perturbations,' Class. Quant. Grav. 37 , 075001 (2020), arXiv:1908.09095 [gr-qc].\n- [76] Manuela Campanelli and Carlos O. Lousto, 'Second order gauge invariant gravitational perturbations of a Kerr  \nblack hole,' Phys. Rev. D 59 , 124022 (1999), arXiv:grqc/9811019.  \n- [77] B. F. Whiting and L. R. Price, 'Metric reconstruction from Weyl scalars,' Class. Quant. Grav. 22 , S589-S604 (2005).\n- [78] John M. Stewart, 'Hertz-Bromwich-Debye-WhittakerPenrose Potentials in General Relativity,' Proc. Roy. Soc. Lond. A 367 , 527-538 (1979).\n- [79] Carlos O. Lousto and Bernard F. Whiting, 'Reconstruction of black hole metric perturbations from Weyl curvature,' Phys. Rev. D 66 , 024026 (2002), arXiv:grqc/0203061.\n- [80] Maarten van de Meent and Abhay G. Shah, 'Metric perturbations produced by eccentric equatorial orbits around a Kerr black hole,' Phys. Rev. D 92 , 064025 (2015), arXiv:1506.04755 [gr-qc].\n- [81] Amos Ori, 'Reconstruction of inhomogeneous metric perturbations and electromagnetic four potential in Kerr space-time,' Phys. Rev. D 67 , 124010 (2003), arXiv:grqc/0207045.\n- [82] David Kofroˇn, 'Point particles and Appell's solutions on the axis of a Kerr black hole for an arbitrary spin in terms of the Debye potentials,' Phys. Rev. D 101 , 064027 (2020), arXiv:2003.10969 [gr-qc].\n- [83] Steffen Aksteiner and Thomas Backdahl, 'Symmetries of linearized gravity from adjoint operators,' J. Math. Phys. 60 , 082501 (2019), arXiv:1609.04584 [gr-qc].\n- [84] E Deadman and J M Stewart, 'Linearized perturbations  \n- of the kerr spacetime and outer boundary conditions in numerical relativity,' Classical and Quantum Gravity 28 , 015003 (2010).\n- [85] S. A. Teukolsky and W. H. Press, 'Perturbations of a rotating black hole. III - Interaction of the hole with gravitational and electromagnet ic radiation,' Astrophys. J. 193 , 443-461 (1974).\n- [86] Edward W Leaver, 'Solutions to a generalized spheroidal wave equation: Teukolsky's equations in general relativity, and the two-center problem in molecular quantum mechanics,' Journal of mathematical physics 27 , 12381265 (1986).\n- [87] Tim Johannsen, 'Photon Rings around Kerr and Kerrlike Black Holes,' Astrophys. J. 777 , 170 (2013), arXiv:1501.02814 [astro-ph.HE].\n- [88] Pierre Christian and Chi-Kwan Chan, 'FANTASY: UserFriendly Symplectic Geodesic Integrator for Arbitrary Metrics with Automatic Differentiation,' Astrophys. J. 909 , 67 (2021), arXiv:2010.02237 [gr-qc].\n- [89] Thomas Madler and Jeffrey Winicour, 'Kerr Black Holes and Nonlinear Radiation Memory,' Class. Quant. Grav. 36 , 095009 (2019), arXiv:1811.04711 [gr-qc].\n- [90] Jos'e A. Font, Jos'e Maria Ib'a˜nez, and Philippos Papadopoulos, 'Non-axisymmetric relativistic Bondi-Hoyle accretion on to a Kerr black hole,' Monthly Notices of the Royal Astronomical Society 305 , 920936 (1999), https://academic.oup.com/mnras/articlepdf/305/4/920/18634678/305-4-920.pdf.", '1. Newman-Penrose formalism': 'For completeness, we review the Newman-Penrose (NP) formalism in this section, adopting the convention from Refs. [68-70]. Note the overall sign difference in some NP quantities compared to the common convention, as seen, for example, in Ref. [16]. This ensures that the background values remain consistent with cases where the metric signature is (+ , -, -, -). In the NP formalism, four null vectors e a ( µ ) are defined as  \ne a ( µ ) = { e a (0) , e a (1) , e a (2) , e a (3) } ≡ { l a , n a , m a , ¯ m a } (A1)  \nwhere l a and n a are real null vectors, and m a and ¯ m a are complex null vectors. They satisfy  \nl a n a = -1 and m a ¯ m a = 1 . (A2)  \nThe metric can be obtained via  \ng ab = -2 l ( a n b ) +2 m ( a ¯ m b ) . (A3)  \nThe directional covariant derivatives are defined as  \nˆ D ≡ l a ∇ a , ˆ ∆ ≡ n a ∇ a , ˆ δ ≡ m a ∇ a , ˆ ¯ δ ≡ ¯ m a ∇ a . (A4)  \nThe Ricci rotation coefficients are given by  \nγ ( a )( b )( c ) ≡ g µλ e µ ( a ) e ν ( c ) ∇ ν e λ ( b ) , (A5)  \ncan be written as [71-73]  \nL † = -L-∇ a e a ( µ ) (A8)  \nTherefore, the formal adjoint directional derivative operators are  \nˆ D † = -( ˆ D + ε + ¯ ε -ρ -¯ ρ ) , ˆ ∆ † = -( ˆ ∆ -γ -¯ γ + µ + ¯ µ ) , ˆ δ † = -( ˆ δ + β -¯ α -τ + ¯ ϖ ) , ˆ ¯ δ † = -( ˆ ¯ δ + ¯ β -α -¯ τ + ϖ ) . (A9)  \nIn the Kerr background, we use the Kinnersley tetrad as follows:  \ne µ (1) = l µ = 1 ∆ BL ( r 2 + a 2 , ∆ BL , 0 , a ) , e µ (2) = n µ = 1 2Σ BL ( r 2 + a 2 , -∆ BL , 0 , a ) , e µ (3) = m µ = -¯ ρ √ 2 ( ia sin θ, 0 , 1 , i sin θ ) , e µ (4) = ¯ m µ = ρ √ 2 ( ia sin θ, 0 , -1 , i sin θ ) , (A10)  \nwhere ρ = -1 / ( r -ia cos θ ) and the overbar denotes complex conjugation. ∆ BL and Σ BL are defined as  \n∆ BL ≡ r 2 -2 Mr + a 2 , Σ BL ≡ r 2 + a 2 cos 2 θ . (A11)  \nThe Ricci rotation coefficients can be written explicitly as  \nκ = λ = ν = σ = ϵ = 0 , (A12)  \nand  \nρ = -1 / ( r -ia cos θ ) , β = -¯ ρ cot θ/ (2 √ 2) , ϖ = iaρ 2 sin θ/ √ 2 , τ = -iaρ ¯ ρ sin θ/ √ 2 , µ = ρ 2 ¯ ρ ∆ / 2 , γ = µ + ρ ¯ ρ ( r -M ) / 2 , α = ϖ -¯ β . (A13)  \nThe Wely scalars are defined as  \nΨ 0 ≡ C (1)(3)(1)(3) = C µνρσ l µ m ν l ρ m σ , Ψ 1 ≡ C (1)(2)(1)(3) = C µνρσ l µ n ν l ρ m σ , Ψ 2 ≡ C (1)(3)(4)(2) = C µνρσ l µ m ν ¯ m ρ n σ , Ψ 3 ≡ C (1)(2)(4)(2) = C µνρσ l µ n ν ¯ m ρ n σ , Ψ 4 ≡ C (2)(4)(2)(4) = C µνρσ n µ ¯ m ν n ρ ¯ m σ , (A14)  \nwhere C µνρσ is the Weyl tensor. In the Kerr background, we have Ψ 0 = Ψ 1 = Ψ 3 = Ψ 4 = 0 and  \nΨ 2 = ρ 3 M. (A15)  \nOther useful equations derived from the Bianchi identities are  \nˆ D Ψ 2 = 3 ρ Ψ 2 , ˆ ∆Ψ 2 = -3 µ Ψ 2 , ˆ δ Ψ 2 = 3 τ Ψ 2 , ˆ ¯ δ Ψ 2 = -3 ϖ Ψ 2 . (A16)  \nwhich can be explicitly written as  \nκ ≡ -γ (3)(1)(1) , ϖ ≡ -γ (2)(4)(1) , ϵ ≡ -1 2 ( γ (2)(1)(1) + γ (3)(4)(1) ) , τ ≡ -γ (3)(1)(2) , ν ≡ -γ (2)(4)(2) , γ ≡ -1 2 ( γ (2)(1)(2) + γ (3)(4)(2) ) , σ ≡ -γ (3)(1)(3) , µ ≡ -γ (2)(4)(3) , β ≡ -1 2 ( γ (2)(1)(3) + γ (3)(4)(3) ) , ρ ≡ -γ (3)(1)(4) , λ ≡ -γ (2)(4)(4) , α ≡ -1 2 ( γ (2)(1)(4) + γ (3)(4)(4) ) . (A6)  \nThe formal adjoint L † of an operator  \nL = e a ( µ ) ∇ a (A7)', '2. Metric Reconstruction': "Cohen, Chrzanowski, and Kegeles (CCK) formulated a method for reconstructing vacuum metric perturbations within a radiation gauge [23-25]. Here, we briefly review the procedure of CCK reconstruction, which is detailed in Ref. [69, 74-78]. Many of our calculations heavily depend on the Newman-Penrose (NP) formalism, adhering to the conventions outlined in Appendix A 1.  \nWe begin by introducing four operators: ˆ E , ˆ O , ˆ S , and ˆ T , each with their respective formal adjoints:  \n- · The linearized Einstein operator ˆ E , which is defined as  \nˆ E ab ( h ) = 1 2 [ -∇ c ∇ c h ab -∇ a ∇ b h c c +2 ∇ c ∇ ( a h b ) c + g ab ( ∇ c ∇ c h d d -∇ c ∇ d h cd )] . (A17)  \nThe operator ˆ E is self-adjoint, as indicated by ˆ E = ˆ E † , and it maps a metric perturbation to its corresponding linearized Einstein tensor.  \n- · The Teukolsky operator ˆ O , which can be written as [16]  \nˆ O 0 = ( ˆ D -3 ϵ +¯ ϵ -4 ρ -¯ ρ )( ˆ ∆ -4 γ + µ ) -( ˆ δ + ¯ ϖ -¯ α -3 β -4 τ )( ˆ ¯ δ + ϖ -4 α ) -3Ψ 2 , ˆ O 4 = ( ˆ ∆+3 γ -¯ γ +4 µ + ¯ µ )( ˆ D +4 ϵ -ρ ) -( ˆ δ -¯ τ + ¯ β +3 α +4 ϖ )( ˆ δ -τ +4 β ) -3Ψ 2 . (A18)  \nThen the Teukolsky equation with spin ± 2 can be expressed as  \nˆ O 0 Ψ (1) 0 = 8 πT 0 , ˆ O 4 Ψ (1) 4 = 8 πT 4 . (A19)  \nTheir separable form, Teukolsky master equation, can be written as  \nˆ O s ψ s = 8 π Σ T s , (A20)  \nwhere the spin s is designated as either ± 2. For s = -2,  \nψ -2 = ρ -4 Ψ (1) 4 , ˆ O -2 = 2Σ ρ -4 ˆ O 4 ρ 4 , T -2 = 2 ρ -4 T 4 ; (A21)  \nconversely, when s = +2,  \nψ +2 = Ψ (1) 0 , ˆ O +2 = 2Σ ˆ O 0 , T +2 = 2 T 0 . (A22)  \n- · The operator ˆ S , which acts on a energy-momentum tensor T ab to yield the source T 4 of the Teukolsky equation, (Eq. 2.15 in Ref. [16])  \nˆ S ab = ( ˆ ∆+3 γ -¯ γ +4 µ + ¯ µ )[ ( ˆ δ -2¯ τ +2 α ) n ( a ¯ m b ) -( ˆ ∆+2 γ -2¯ γ + ¯ µ ) ¯ m a ¯ m b ] + ( ˆ δ -¯ τ + ¯ β +3 α +4 ϖ )[( ˆ ∆+2 γ +2¯ µ ) n ( a ¯ m b ) -( ˆ δ -¯ τ +2 ¯ β +2 α ) n a n b ] , (A23)  \n- · The operator ˆ T , which acts a metric perturbation h ab to the perturbed Weyl scalar Ψ (1) 4 .  \nˆ T ab = -1 2 {( ˆ δ -¯ τ +3 α + ¯ β )( ˆ δ -¯ τ +2 α +2 ¯ β ) n a n b + ( ˆ ∆+¯ µ +3 γ -¯ γ )( ˆ ∆+¯ µ +2 γ -2¯ γ ) ¯ m a ¯ m b -[( ˆ ∆+¯ µ +3 γ -¯ γ ) ( ˆ δ -2¯ τ +2 α ) + ( ˆ δ -¯ τ +3 α + ¯ β )( ˆ ∆+2¯ µ +2 γ )] n ( a ¯ m b ) } (A24)  \nAs pointed out by Wald [71], one can find the following operator identity  \nˆ O 4 ˆ T = ˆ S ˆ E , (A25)  \nh ab ˆ E ˆ E † T ab ˆ S ˆ S † T 4 ˆ O 4 ˆ O † 4 Ψ (1) 4 ˆ T † ˆ T  \nFIG. 1: Relation of operators for s = -2.  \nand its formal adjoint, ˆ T † ˆ O † 4 = ˆ E ˆ S † , is also valid. Then, assuming a complex potential Φ, also known as Herz potential, satisfies the adjoint Teukolsky equation ˆ O † 4 Φ = 0, we can consequently have  \nˆ E ab ( h ) = 2 Re( ˆ T † ab ˆ O † 4 Φ) = 0 (A26)  \nTherefore, under the assumption, h ab provides a solution to the linearized Einstein equation, and can be expressed as  \nh ab = 2Re ( ˆ S † ab Φ ) . (A27)  \nChrzanowski introduced two radiation gauges: the ingoing radiation gauge (IRG) and the outgoing radiation gauge (ORG) [23]. In this paper, we choose ORG, which can be written as  \nh ab n b = 0 = g ab h ab . (A28)  \nUnder this gauge condition, the metric satisfies  \nh nn = 0 = h ln = 0 = h nm = 0 = h n ¯ m = 0 = h m ¯ m . (A29)  \nThe Hertz potential Φ ORG fulfills the adjoint Teukolsky equation O † 4 Φ ORG = 0, where  \nO † 4 = ( ˆ D -3 ϵ +¯ ϵ -¯ ρ )( ˆ ∆ -4 γ -3 µ ) -( ˆ δ -3 β -¯ α + ¯ ϖ )( ˆ ¯ δ -4 α -3 ϖ ) -3Ψ 2 . (A30)  \nOne can verify O † 4 Φ ORG = 0 is equivalent to the homogeneous Teukolsky equation for a field with s = +2, that is,  \nˆ O † 4 Φ ORG = ρ -4 O 0 ( ρ 4 Φ ORG ) = 0 , (A31)  \nand the metric component can be written as [72]  \nh ab = { n ( a ¯ m b ) [( ˆ ¯ δ + ¯ β -3 α + ¯ τ + ϖ )( ˆ ∆ -4 γ -3 µ ) + ( ˆ ∆ -3 γ -¯ γ + µ -¯ µ )( ˆ ¯ δ -4 α -3 ϖ )] -n a n b ( ˆ ¯ δ -¯ β -3 α + ϖ )( ˆ ¯ δ -4 α -3 ϖ ) -¯ m a ¯ m b ( ˆ ∆ -3 γ + ¯ γ + µ )( ˆ ∆ -4 γ -3 µ )} Φ ORG +c.c. , (A32)  \nwhere 'c.c.' stands for complex conjugate part of the whole object. If we define ˚ Φ ORG = ρ 4 Φ ORG , then the non-zero metric component can be written as [79, 80]  \nh ll ≡ h µν l µ l ν = ˆ H † ll ˚ Φ ORG +c.c. , h lm ≡ h µν l µ m ν = ˆ H † lm ˚ Φ ORG , h l ¯ m ≡ ¯ h lm , h mm ≡ h µν m µ m ν = ˆ H † mm ˚ Φ ORG , (A33)  \nwhere  \nˆ H † ll = -ρ -4 ( ˆ ¯ δ -3 α -¯ β +5 ϖ )( ˆ ¯ δ -4 α + ϖ ) , ˆ H † lm = -ρ -4 2 {( ˆ ¯ δ -3 α + ¯ β +5 ϖ + ¯ τ )( ˆ ∆+ µ -4 γ ) + ( ˆ ∆+5 µ -¯ µ -3 γ -¯ γ )( ˆ ¯ δ -4 α + ϖ )} , ˆ H † mm = -ρ -4 ( ˆ ∆+5 µ -3 γ + ¯ γ )( ˆ ∆+ µ -4 γ ) . (A34)  \nWhile ˚ Φ ORG does indeed satisfy the vacuum Teukolsky equation with s = +2, the field does not correspond with either the Weyl scalar Ψ (1) 0 or Ψ (1) 4 originating from the reconstructed metric h ab . From Fig. 1, it is clear that Ψ (1) 4 can be obtained via a fourth-order ordinary differential equation (ODE) along ingoing null rays [69, 80, 81]  \nΨ (1) 4 = ˆ T ab h ab . (A35)  \nWe can reduce this circularity condition to what is known as a radial inversion relation [69, 70, 82-84]  \nΨ (1) 4 = 1 2 ( ˆ ∆+3 γ -¯ γ )( ˆ ∆+2 γ -2¯ γ )( ˆ ∆+ γ -3¯ γ )( ˆ ∆ -4¯ γ ) ¯ Φ ORG , (A36)  \nwhich can be solved in vacuum algebraically by applying the separation of variable in next section.", 'a. QNM to Hertz potential': "In Boyer-Lindquist (BL) coordinates, the Teukolsky operator ˆ O s with spin s in Eq. (A20) can be written as [16]  \nˆ O s ψ s = [ ( r 2 + a 2 ) 2 ∆ BL -a 2 sin 2 θ ] ∂ 2 ψ s ∂t 2 + 4 Mar ∆ BL ∂ 2 ψ s ∂t∂ϕ + [ a 2 ∆ -1 sin 2 θ ] ∂ 2 ψ s ∂ϕ 2 -∆ -s BL ∂ ∂r ( ∆ s +1 BL ∂ψ s ∂r ) -1 sin θ ∂ ∂θ ( sin θ ∂ψ s ∂θ ) -2 s [ a ( r -M ) ∆ BL + i cos θ sin 2 θ ] ∂ψ s ∂ϕ -2 s [ M ( r 2 -a 2 ) ∆ BL -r -ia cos θ ] ∂ψ s ∂t + ( s 2 cot 2 θ -s ) ψ s (A37)  \nwhere ψ s is a field with spin s . In the vacuum case T s = 0, we can apply the separation of variables by letting  \nψ s ( t, r, θ, ϕ ) = 1 √ 2 π ∑ l,m,n s Z ℓmn s S ℓmn ( θ ) s R ℓmn ( r ) e -iω mn t + imϕ , (A38)  \nWith z = cos θ , s S ℓmn ( θ ) satisfy  \n∂ z [( 1 -z 2 ) ∂ z [ s S ℓmn ( θ )] ] + [ ( cz ) 2 -2 csz + s + A ℓmn s ( c ) -( m + sz ) 2 1 -x 2 ] s S ℓmn ( θ ) = 0 (A39)  \nwhere c = aω and s S ℓmn ( θ ) e imϕ is the spin-weighted spheroidal harmonics. To normalize the spin-weighted spheroidal harmonics, we apply the following normalization condition  \n∫ π 0 s S ℓmn ( θ ) 2 sin θ d θ = 1 . (A40)  \nThe radial function R ( r ) is required to fulfill  \n∆ -s BL d dr [ ∆ s +1 BL d s R ℓmn ( r ) d r ] + [ K 2 -2 is ( r -M ) K ∆ BL +4 isωr -λ ℓmn ] s R ℓmn ( r ) = 0 . (A41)  \nwhere  \nK mn ≡ ( r 2 + a 2 ) ω mn -am, s λ ℓmn ≡ s A ℓmn + a 2 ω 2 mn -2 amω mn . (A42)  \nHere s A ℓmn is the angular separation constant.  \nIn the following, we will revisit the procedure for obtaining the Hertz potential corresponding to specific Weyl scalar Ψ (1) 0 or Ψ (1) 4 , which is discussed in Ref. [80, 81]. According to the previously mentioned condition, ˆ O +2 ˚ Φ ORG = 0, we can expand ˚ Φ ORG in modes  \n˚ Φ ORG = 1 √ 2 π ∑ ℓ,m,n C ℓmn +2 R ℓmn ( r ) +2 S ℓmn ( θ ) e imϕ -iω mn t , (A43)  \nwhere C ℓmn are constant coefficients. In addition, the Hertz potential ˚ Φ ORG also satisfies Eq. (A36), which can be simplified using the following relationship  \n( ˆ ∆ -m ¯ γ + nγ ) f = ∆ ( n -m ) / 2 BL ¯ ρ -m ρ n ˆ ∆ [ ∆ ( m -n ) / 2 BL ¯ ρ m ρ -n f ] . (A44)  \nIn this equation, f = f ( t, r, θ, ϕ ) can represent any function. Then Eq. (A36) can be written as [70, 79-81]  \n1 32 ∆ 2 BL ( ˆ D † 0 ) 4 [ ∆ 2 BL ¯ ˚ Φ ORG ] = ρ -4 Ψ (1) 4 . (A45)  \nHere, the operator ˆ D † 0 is defined as:  \nˆ D † 0 = ∂ r -( r 2 + a 2 ) ∂ t + a∂ ϕ ∆ BL . (A46)  \nNote that Eq. (A45) can be solved easily on the mode basis given by Eq. (A43), using the operator  \nˆ D † mn = ∂ r + i K mn ∆ BL = ∂ r + i ω mn ( r 2 + a 2 ) -ma ∆ BL , (A47)  \nwhich act on radial modes. The coefficient can be fixed by considering the asymptotic limit approaching infinity and the horizon. This analysis was previously carried out in Ref. [80, 81] for real ω . Due to our choice of gauge (ORG), the results are similar to those in Section III.C of Ref. [80]. Here, we briefly review and extend the procedure to complex ω .  \n1 32 ∆ 2 BL ( ˆ D † 0 ) 4 [ ∆ 2 BL ¯ ˚ Φ ORG ] = ∑ ℓ,m,n ¯ C ℓmn 32 √ 2 π ∆ 2 BL ( ∂ r -i ¯ K mn ∆ BL ) 4 [ ∆ 2 BL +2 ¯ R ℓmn ( r ) +2 S ℓmn ( θ ) e -imϕ + i ¯ ω mn t ] (A48)  \nBy relabelling ( m,n ) to ( -m, -n ) and employing the identities s S ℓmn ( θ ) = ( -1) s + m -s S l -m -n ( θ ), s ¯ R ℓmn = s R l -m -n , ¯ K mn = -K -m -ns and ¯ ω mn = -ω -m -n , the equation transforms into  \n1 32 ∆ 2 BL ( ˆ D † 0 ) 4 [ ∆ 2 BL ¯ ˚ Φ ORG ] = ∑ ℓ,m,n ( -1) m ¯ C ℓ -m -n 32 √ 2 π ∆ 2 BL ( ˆ D † mn ) 4 [ ∆ 2 BL +2 ¯ R ℓmn ( r ) -2 S ℓmn ( θ ) e imϕ -iω mn t ] (A49)  \nWhen considering Eq. (A38), we obtain  \n-2 Z ℓmn -2 R ℓmn ( r ) = ( -1) m 32 ¯ C ℓ -m -n ∆ 2 BL ( ˆ D † mn ) 4 [ ∆ 2 BL +2 R ℓmn ( r ) ] . (A50)  \nThe asymptotic behavior of s R ℓmn is as follows [20, 85]:  \nlim r ∗ →∞ s R ℓmn ( r ∗ ) ∼      e -iωr ∗ r : ingoing wave e iωr ∗ r 2 s +1 : outgoing wave. (A51)  \nlim r ∗ →-∞ s R ℓmn ( r ∗ ) ∼ { e ik mn r ∗ : ingoing (out of BH) ∆ -s BL e -ik mn r ∗ : outgoing (into BH) (A52)  \nwhere k mn = ω mn -ma 2 Mr + , and r ∗ is the tortoise coordinate which is given by  \nr ∗ = r + 2 M r + -r -[ r + ln ( r -r + 2 M ) -r -ln ( r -r -2 M )] (A53)  \nIf f ( r ) is a solution to the radial Teukolsky equation Eq. (A41) with spin s , then ∆ s BL ¯ f ( r ) is a solution to the radial Teukolsky equation with spin -s . For generic modes with ω mn = 0, the asymptotic behavior suggests that this  \n̸  \nrelation swaps physical (retarded) boundary conditions with unphysical (advanced) ones. We denoted the relationship as  \n-s R # ℓmn ≡ ∆ s BL s ¯ R ℓmn . (A54)  \nWith this in mind, we have  \n-2 Z ℓmn -2 R ℓmn ( r ) = ( -1) m 32 ¯ C ℓ -m -n ∆ 2 BL ( ˆ D † mn ) 4 [ -2 R # ℓmn ] . (A55)  \nwhen we insert the asymptotic behavior into Eq. (A55), we can find  \nC ℓmn = 2 ω 4 mn ( -1) ℓ + m -2 Z ℓmn (A56)  \nHere, we used the following relationship  \n-2 ¯ Z lmn = ( -1) l -2 Z l -m -n . (A57)  \nConsequently, once we have the solution set { ω, Ψ (1) 4 } for the Teukolsky equation with s = -2, we are able to reconstruct the first-order metric perturbation h ab . After reconstructing the metric, we verify its validity by substituting it into the linearized Einstein equation ˆ E ab ( h ) = 0.  \nWecalculate the quasinormal mode using the qnm package [36], which employs Leaver's method to the radial equation and the Cook-Zalutskiy spectral approach to the angular sector [20, 86]. Furthermore, we use the pseudospectral method to get solutions to the radial Teukolsky equation. Here we briefly illustrate our steps. We convert the asymptotic behavior into coordinate r [20]  \nlim r →∞ s R ℓmn ( r ) ∼ r -1 -2 s +2 iω mn M e iω mn r (A58)  \nlim r → r + s R ℓmn ( r ) ∼ ( r -r + ) -s -iσ + (A59)  \nwhere σ + = 2 ω mn Mr + -ma r + -r -. Then we can define a new radial function s ˜ R ℓmn ( r ) by  \ns R ℓmn ( r ) = ( r -r + ) -s -iσ + r -1 -s +2 iω mn M + iσ + e iω mn r s ˜ R ℓmn ( r ) , (A60)  \nwhere s ˜ R ℓmn ( r ) is regular at BH horizon and infinity, and we normalize it by applying  \nlim r →∞ s ˜ R ℓmn ( r ) = 1 (A61)  \nWe can expand s ˜ R ℓmn ( r ) as a series of Chebyshev polynomials and get s S ℓmn ( θ ) in terms of spin-weighted spherical functions as in Ref. [20]. As a result, we can use our numerical results and evaluate them for any point ( t, r, θ, ϕ ) of spacetime in the subsequent ray tracing section. To ensure the perturbation metric is flawless, we also verified ˆ E ab ( h ) = 0 numerically.", '3. Radiated energy': 'In this section, we estimate the radiated energy of a single mode after truncation. The energy flux is given by [34]  \nd E d t = lim r →∞ [ r 2 16 π ∫ Ω ∣ ∣ ∣ ∣ ∫ t + ∞ Ψ (1) 4 d ˜ t ∣ ∣ ∣ ∣ 2 dΩ ] , (A62)  \nThen the total energy emitted can be calculated as  \nE tot = ∫ ∞ t 0 d E d t d t (A63)  \nwhere t 0 = r ∗ ( r o ). This total energy satisfies the condition E tot /M ⪅ 3% [34]. Assuming that only the ℓ = m = 2 mode dominates the radiated energy, this relationship leads to  \n| s Z | ⪅ 0 . 24 . (A64)  \nIt is important to note that this approach differs slightly from truncating the metric using H ( t -r ∗ ). However, it should provide a sufficiently accurate estimation for our purposes.', '1. Local oberver basis': 'The observation of photons at any given point is dependent on the observer. For the case of a pure Kerr black hole, zero angular momentum observers (ZAMOs) are typically employed. In this section, we extend them to dynamically perturbed Kerr spacetime.  \nWithin the Boyer-Lindquist coordinate system, a ZAMO frame can be chosen and expressed in terms of the coordinate basis { ∂ t , ∂ r , ∂ θ , ∂ ϕ } [26, 27, 30, 87] as  \n⃗z µ ( t ) = g ϕϕ ∂ t -g ϕt ∂ ϕ √ g ϕϕ ( g 2 ϕt -g ϕϕ g tt ) , ⃗z µ ( r ) = ∂ r √ g rr , ⃗z µ ( θ ) = ∂ θ √ g θθ , ⃗z µ ( ϕ ) = ∂ ϕ √ g ϕϕ . (B1)  \nThe observer basis ⃗z ν ( σ ) has a Minkowski normalization  \ng (0) µν ⃗z µ ( ρ ) ⃗z ν ( σ ) = η ( ρ )( σ ) , (B2)  \nwith η ( ρ )( σ ) = diag {-1 , 1 , 1 , 1 } . Unlike the unperturbed case, Eq. (B2) no longer holds when the unperturbed metric g (0) µν is replaced by the perturbed metric g µν . To address this issue, we employ a Gram-Schmidt orthogonalization process, using the chosen ZAMO tetrad ⃗v ρ = ⃗z µ ( ρ ) as the initial set of vectors. The Gram-Schmidt process can be expressed as  \n⃗u 0 = ⃗v 0 , ⃗u 1 = ⃗v 1 -proj ⃗u 0 ( ⃗v 1 ) , ⃗u 2 = ⃗v 2 -proj ⃗u 0 ( ⃗v 2 ) -proj ⃗u 1 ( ⃗v 2 ) , ⃗u 3 = ⃗v 3 -proj ⃗u 0 ( ⃗v 3 ) -proj ⃗u 1 ( ⃗v 3 ) -proj ⃗u 2 ( ⃗v 3 ) , (B3)  \nwhere the projection operator is defined as  \nproj ⃗u ⃗v = ⟨ ⃗u, ⃗v ⟩ ⟨ ⃗u, ⃗u ⟩ ⃗u , (B4)  \nwith the inner product evaluated using the metric g µν . The resulting orthonormalized tetrad, denoted as  \n˜ ⃗z µ ( ρ ) = ⃗u ρ √ |⟨ ⃗u ρ , ⃗u ρ ⟩| (B5)  \nsatisfies the following condition  \ng µν ˜ ⃗z µ ( ρ ) ˜ ⃗z ν ( σ ) = η ( ρ )( σ ) . (B6)', '2. Impact Parameters': "In the image plane of the observer defined above, each photon is assigned Cartesian coordinates ( ˜ X, ˜ Y ), which represent its impact parameters. These coordinates are proportional to the photon's observation angles (˜ α, ˜ β ). Here we follow a convention similar to Ref. [30], and define them as:  \n˜ X ≡ -r ˜ β, ˜ Y ≡ r ˜ α. (B7)  \nConsidering the geometry of photon detection, we can express the photon's 4-momentum in the observer's reference frame as follows (refer to Fig. 1 in Ref. [30]):  \np ( ϕ ) = | ⃗ P | sin ˜ β cos ˜ α, p ( θ ) = | ⃗ P | sin ˜ α, p ( r ) = | ⃗ P | cos ˜ β cos ˜ α, (B8)  \nwhere | ⃗ P | = p ( t ) . From this, we can derive the photon's momentum using  \np µ = ˜ ⃗z µ ( σ ) p ( σ ) . (B9)", '3. Kerr-Schild coordinates': 'To avoid the numerical problem when a photon gets too close to the pole in the BL coordinate system, we perform ray tracing in the Cartesian Kerr-Schild (KS) coordinate system { t KS , x, y, z } . The coordinate transformations, which can be found in Ref. [88-90], are also reviewed in the following.  \nThe Kerr metric in Kerr-Schild coordinates { t KS , x, y, z } can be expressed as  \ng (0) µν = η µν + fk µ k ν , f = 2 Mr 3 r 4 + a 2 z 2 , k µ = ( 1 , rx + ay r 2 + a 2 , ry -ax r 2 + a 2 , z r ) , (B10)  \nwhere r is the same as the radial coordinate in the BL coordinates. The radial coordinate r can be defined in terms of the Cartesian KS coordinates as  \nx 2 + y 2 r 2 + a 2 + z 2 r 2 = 1 . (B11)  \nThen the coordinates { r, θ, ϕ KS } can be written explicitly in terms of Cartesian KS as  \nr =   R 2 -a 2 + √ ( R 2 -a 2 ) 2 +4 a 2 z 2 2   1 / 2 , θ = arccos z r , ϕ KS = tan -1 ( -ay -rx, ax -ry ) + π . (B12)  \nwhere  \nR = √ x 2 + y 2 + z 2 . (B13)  \nThe function tan -1 ( x, y ), also known as the two-argument arctangent function, calculates the arctangent of y/x while considering the quadrant in which the point ( x, y ) lies.  \nThe coordinate transformation from BL coordinates can be written as  \nt KS = t + ∫ 2 Mr ∆ BL dr , x = √ r 2 + a 2 sin θ cos [ ϕ KS +arctan( a/r )] , y = √ r 2 + a 2 sin θ sin [ ϕ KS +arctan( a/r )] , z = r cos θ , (B14)  \nwhere  \nϕ KS = ϕ + a ∫ dr ∆ BL . (B15)  \nFollowing Ref. [88], we fix the integration constant is fixed by letting t KS = t and ϕ KS = ϕ at the initial radius r o , so we have  \nt KS = t + M √ M 2 -a 2 [ r + ln ( r -r + r o -r + ) -r -ln ( r -r -r o -r -)] , ϕ KS = ϕ + a 2 √ M 2 -a 2 [ ln ( r -r + r o -r + ) -ln ( r -r -r o -r -)] . (B16)'}