Frequency- and polarization-dependent lensing of gravitational waves in strong gravitational fields

The propagation of gravitational waves can be described in terms of null geodesics by using the geometrical optics approximation. However, at large but finite frequencies the propagation is affected by the spin-orbit coupling corrections to geometrical optics, known as the gravitational spin Hall effect. Consequently, gravitational waves follow slightly different frequency- and polarization-dependent trajectories, leading to dispersive and birefringent phenomena. We study the potential for detecting the gravitational spin Hall effect in hierarchical triple black hole systems, consisting of an emitting binary orbiting a more massive body acting as a gravitational lens. We calculate the difference in time of arrival with respect to the geodesic propagation and find that it follows a simple power-law dependence on frequency with a fixed exponent. We calculate the gravitational spin Hall-corrected waveform and its mismatch with respect to the original waveform. The waveform carries a measurable imprint of the strong gravitational field if the source, lens, and observer are sufficiently aligned, or for generic observers if the source is close enough to the lens. We present constraints on dispersive time delays from GWTC-3, translated from limits on Lorentz invariance violation. Finally, we address the sensitivity of current and future ground detectors to dispersive lensing. Our results demonstrate that the gravitational spin Hall effect can be detected, providing a novel probe of general relativity and the environments of compact binary systems.


I. INTRODUCTION
The first detection of gravitational waves (GWs) -GW150914 -by the Advanced LIGO observatory marked the beginning of the new era of gravitational-wave astronomy [1,2].GWs carry information about their source, but also imprints of the spacetime on which they travel.Observable sources of GWs emit over a wide range of frequencies [3].As an example, the aforementioned GW150914 was detected from ∼ 35 to 250 Hz.Its wavelength and that of any signal detectable by LIGO-Virgo-Kagra (LVK), remains orders of magnitude larger than that of the longest electromagnetic (EM) signal capable of crossing the atmosphere (∼ 10 MHz).Therefore, GWs have the potential to detect novel propagation effects at low frequencies, particularly when their wavelength approaches characteristic lengths of physical systems -e.g. the Schwarzschild radius of a black hole or other gravitational lens.In these cases, the propagation of GWs might deviate slightly from the standard predictions of geometrical optics (GO) [4][5][6].
The GSHE is described by a set of effective ray equations that represent the propagation of a gravitational wave packet energy centroid up to first order in wavelength, derived as a higher-order GO approximation using a Wentzel-Kramers-Brillouin (WKB) ansatz.Within this formalism, the wave packets undergo frequency-and polarization-dependent deviations from the GO trajectory, which can be viewed as a manifestation of the spinorbit coupling via the Berry curvature.Moreover, the deviations are described by the same effective ray equations for both EM and linearized gravitational fields [22,24]. 1 In this paper, spin-orbit coupling (see, e.g., Ref. [7]) refers to the dynamics of wave packets with internal structure, where the spin represents the internal degree of freedom of the wave packet (i.e., polarization), while the orbital part refers to the motion of the wave packet as a whole.Thus, this should not be confused with spin-orbit couplings arising in the dynamics of black hole binaries during the coalescence process [8].
GWs offer the best chance to probe the GSHE.The GSHE emerges as a first-order perturbation in the ratio between wavelength and the background gravitational field length scale -the Schwarzschild radius R s .While the present day GW terrestrial observatories have a lower limit at ∼ 10 Hz [45], or equivalently wavelength of ∼ 10 7 m, radio telescopes such as the Event Horizon telescopes observe at ∼ 1.3×10 −3 m [46], i.e. at wavelengths orders of magnitude lower than the GW interferometers.Therefore, there is little hope of finding observable astrophysical systems where the EM radiation wavelength is comparable to the gravitational field length scale.
Another reason to search for GSHE using GWs is that sources may inhabit high-curvature environments.In addition to isolated evolution of massive binary stars, GW emitting binaries may form by dynamical encounters in a dense environment, such as a globular cluster [47][48][49] or an AGN [50][51][52].For a review of hierarchical black hole (BH) formation channels, see Ref. [53].In the active galactic nucleus (AGN) scenario, compact objects accumulate in the disk around a supermassive black hole [54].Interactions with the disk would subsequently drive them towards migration traps, stable orbits where gas torques change direction [55].Migration traps could be as close as ≲ 10 Schwarzschild radii of the supermassive black hole [56].Such a "last migration trap" may contribute up to ∼ 1% of GW events detectable by LVK.This opens up the possibility of detecting strong field effects in GW propagation in hierarchical triple systems, wherein the emitting BH binary is sufficiently close to or orbiting around a massive third companion BH.The GSHE may be detectable in these systems, in addition to multiple images of the merger caused by the BH [57][58][59].
Interest in the AGNs-GW connection boomed after LIGO-Virgo's detection of GW190521 [60,61], a binary whose primary component's mass is in the pair instability gap [62].Such a massive BH could not have formed from stellar evolution, pointing towards a likely dynamical origin for the binary.Furthermore, the Zwicky Transient Facility detected an EM flare in AGN J124942.3+344929(redshift of 0.438), 34 days after GW190521 and with consistent sky localization.In this tentative interpretation, the BH binary would be in a migration trap with a semi-major axis of ∼ 350 Schwarzschild radii of the supermassive black hole, and the delay between both events would be the time required for the EM radiation to emerge from the accretion disk of the AGN [63].Although suggestive, evidence for an AGN origin of GW190521 is far from conclusive when considering LIGO-Virgo data [64][65][66][67][68] or the putative EM counterpart [69,70].
The GSHE provides a novel test of the GW source environments, which may help establish their AGN formation channel.An advantage of this test is that it can be performed on individual observations.In contrast, other proposed methods require either LISA-like observatory [71,72] to measure the orbit of the emitting binary around the background black hole [58,[73][74][75][76] or popu-lation studies.The latter being based on binary properties (masses, spin, eccentricity) [77][78][79] or associating GW events with detected AGN flares [70,80].Although measuring the GSHE might be possible for only a small fraction of the GW events originating in AGN disks, its complementarity with other methods would yield valuable insights into BH and GW astrophysics.
The GSHE arises in Einstein's general theory of relativity (GR) [24], but it is also similar to effects emerging in theories beyond GR, and thus needs to be taken into account to correctly interpret tests of gravity with GWs.A nonzero graviton mass leads to a distance and frequency-dependent propagation for all GWs [81].Some alternative theories predict environment-and polarization-dependent GW propagation speeds -the GW birefringence effect [82].This leads to a frequencyindependent time delay between the + and × polarization states that may either interfere in the detector or appear as two copies of the same signal if the time delay is shorter/longer than the signal, respectively.A related effect stems from parity-breaking terms in the effective field theory of GWs.Ref. [83] searched for frequencydependent GW birefringence (between left and right polarized GWs), finding that only the GW190521 observation is compatible with violation of parity.All these beyond-GR effects are related to the GSHE, although in principle distinguishable from it.Establishing a detection of the GSHE in the GW data would represent yet another test of gravity and additional evidence for GR in the strong-field regime.
We demonstrate that the GSHE can be detected in GW sources in a hierarchical triple system, in which a stellar-mass binary is close to a much more massive companion, such as in an AGN.The main observable signature of the GSHE is time delay between the high-and low-frequency components of the waveform, with a correction proportional to ∼ 1/f 2 relative to geodesic propagation.Therefore, the GSHE may appear as an inconsistency between the higher and lower frequency parts of the waveform (e.g., in inspiral-mergerringdown tests of GR).A subdominant GSHE signature is a frequency-dependent birefringence effect -a time delay between the left-and right-polarized components.GSHE-birefringence is further suppressed (∼ 1/f 3 ) and is likely too small to be detectable, except in fine-tuned configurations.A third signature of this scenario is the likely presence of multiple signals due to strong-field lensing by the massive BH.The relative magnification, time delay, and sign of the GSHE correction between these signals should allow for further means to probe the system configuration.
The paper is organized as follows.We begin by describing the GSHE and the numerical calculation of the time of arrival delay in Section II.In Section III, we describe the dependence of the time of arrival delay on frequency, polarization state, and the mutual position of the source and the observer.In Section III, we demonstrate the effect of the GSHE on a GW waveform and its distin-guishability from an uncorrected waveform.Lastly, we discuss our findings in Section IV and conclude in Section V. Our results are also presented in a more compact form in the companion Letter, Ref. [84].
We note that log refers to a logarithm of base 10, x•y = x µ y µ denotes the inner product of 4-vectors and, unless explicitly discussing dimension-full quantities, we set the speed of light, the gravitational constant and the Kerr BH mass M to unity, c = G = M = 1.

II. METHODOLOGY
We assume the existence of a GW emitter -a binary BH merger -in the vicinity of a Kerr BH, with GW ray trajectories passing through the strong-field regime of the background Kerr metric.We then calculate the observer time of arrival of the GSHE trajectories, which depends on frequency and polarization, and compare it to the geodesic time of arrival.In other words, the observer detects that the waveform modes have a frequency-and polarization-dependent time of arrival that deforms the resulting waveform.
We start by reviewing the Kerr metric and GSHE equations in Section II A. We then present our geometric setup in Section II B and numerical integration in Section II C. Finally, we characterize the GSHE time delay quantities in Section II D and discuss our waveform model in Section II E.

A. Gravitational spin Hall equations
We consider the background spacetime of a Kerr black hole with mass M = 1 and spin parameter a, described using Boyer-Lindquist coordinates (t, r, θ, ϕ) [85, p.195].The line element is where We also consider an orthonormal tetrad ) ) ) that satisfies (e a ) µ (e b ) µ = η ab , where η ab is the Minkowski metric.The vectors e a will be used in the definition of the GSHE and for the prescription of initial conditions.On the Kerr background spacetime, we consider GWs represented by small metric perturbations and described by the linearized Einstein field equations.High-frequency GWs can be described using the GO approximation [86,Sec. 35.13], in which case their propagation is determined by the null geodesics of the background spacetime.However, at high but finite frequencies, higher-order corrections to the GO approximation become important.
In this paper, we consider first order in wavelength corrections to the GO approximation, wherein the propagation of GWs is frequency-and polarization-dependent.This is known as the GSHE [24], and the propagation of circularly polarized gravitational wave packets is described by the GSHE equations [24,26] where x µ (τ ) is the worldline of the energy centroid of the wave packet, p µ (τ ) is the average momentum of the wave packet, the spin tensor S αβ describes the angular momentum carried by the wave packet and T α is a timelike vector field with respect to which the energy centroid of the wave packet is defined.We eliminate the ODE for p 0 by enforcing the null momentum condition p • p = 0 along the worldline.For the circularly polarized wave packets that we consider here, the spin tensor is uniquely fixed as where s = ±2, depending on the state of circular polarization.In the context of the high-frequency analysis [24,26] leading to the above equations, the wave frequency f measured by an observer with 4-velocity T α is defined as p • T = −ϵf .The small dimensionless parameter ϵ has the same meaning as in standard high-frequency approximations in general relativity (see, for example, Refs.[87,Sec. 1.5] and [88,Sec. 3.2]), and is meant to keep track of the order of different terms in these expansions.In particular, the GSHE equations were derived in Ref. [24] under the assumption that the wavelength λ is much smaller than the length scale L over which the spacetime varies significantly.Thus, the small expansion parameter is defined as ϵ = O(λ/L), and the GSHE equations provide a reasonable approximation only when ϵ < 1.In the case we are considering here, the lengthscale over which spacetime varies significantly is set by the size of the black hole, so we can take where λ is the wavelength of the wave packet in the rest frame of the source.This can also be expressed in dimension-full quantities as The GSHE equations in Eq. (2.4) depend on the choice of a timelike vector field T α .The role of this vector field has been discussed in detail in Ref. [26], where it has been shown to have physical meaning only at the point of emission and the point of observation of a polarized ray.At these points, T α can be identified with the 4-velocity of the source and observer, respectively, and is responsible for the relativistic Hall effect [89,90].Nevertheless, one has to choose a smooth vector field T α defined everywhere in the region where the GSHE equations are to be integrated.We discuss our choice of T α in the following subsection.

B. Spatial configuration
We consider a static source of GWs close to the BH at x src = (r src , θ src , ϕ src ) with a 4-velocity T α src and a static observer far from the BH at x obs = (r obs , θ obs , ϕ obs ) with a 4-velocity T α obs .The timelike vector field T α appearing in the GSHE equations (2.4) is chosen such that We start with the orthonormal tetrad e a from Eq. (2.3) and perform a spacetime-dependent local Lorentz boost of the orthonormal tetrad such that (e 0 ) α maps to T α src and T α obs at x src and x obs , respectively.We can express the boosted orthonormal tetrad ẽa as where (2.10) The exponential factor ensures a smooth transition between T α src , e α 0 and T α obs .We identify the timelike observer vector field in the GSHE equations (2.4) as T α = (ẽ 0 ) α and further justify the Lorentz boost in Appendix A 2.
For simplicity, we consider a static isotropic emitter of GWs in the vicinity of a massive "lensing" BH that sources the background Kerr metric and a far static observer measuring the waveform (wave packet).The caveat of isotropic emission is relevant as the emission direction of the ϵs-parameterized bundle trajectories must be rotated with respect to the geodesic, ϵ = 0, emission direction by an angle ∼ ϵs.In this work, we do not account for the directional dependence as it is a subdominant effect.Including it would necessitate accounting for it while generating the waveform frequency modes.The Boyer-Lindquist coordinate time t can be related to the static observer's proper time τ as which we derive in Appendix A 1, with g µν being the Kerr metric tensor.Throughout the rest of the paper, we denote the coordinate time as t and the static observer proper time as τ .
A signal with initial momentum p init emitted by the static source with 4-velocity T src has a frequency f src in the source frame.On the other hand, the static observer with 4-velocity T obs will measure the signal frequency as f obs .The observer, therefore, measures the signal redshifted by where p f is the wave packet's momentum when it reaches the observer.This is the common expression for gravitational redshift, which is satisfied up to first order in ϵ.
The ϵ dependence of the gravitational redshift originates from p f and p init , as the initial conditions of a trajectory between two fixed spatial locations depend on ϵ.We will find the ϵ dependence of the gravitational redshift to be negligible.Therefore, since this produces a uniform frequency offset and no new effect, we do not consider this further.Moreover, upon emission, the following relation is enforced, where the last equality follows from Eq. (2.6).An analogue of this condition is then satisfied along the trajectory, as discussed in Ref. [26].We parameterize p init by a unit three-dimensional Cartesian vector k expressed in spherical coordinates where 0 ≤ ψ ≤ π and 0 ≤ ρ < 2π are the polar and azimuthal angle, respectively.The angles ψ and ρ represent the emission direction on the source celestial sphere, and we have that p init = ẽ0 + sin ψ cos ρẽ 1 + sin ψ sin ρẽ 2 + cos ψẽ 3 , (2.14) which satisfies both Eq.(2.13) and the null momentum condition p • p = 0.The initial momentum pointing towards the BH, i.e. with an initial negative radial component, may equally be parameterized with k 2 and k 3 , which can be related to ψ and ρ as (2.16b) We calculate the magnification factor µ defined as the ratio of the source area to the image area [88,91,92].In our case, a trajectory defines a mapping from the celestial sphere of source to the far sphere of radius r = r obs centered at the origin, which we can write as (ψ, ρ) → (θ, ϕ).The magnification µ is where the Jacobian J is defined as The magnification scales a signal propagated along a trajectory by a factor of |µ| and the signal parity is given as the sign of µ, or equivalently the sign of det J. Therefore, as a consequence of the GSHE the magnification is dependent on frequency and polarization.We will explicitly denote this dependence as µ(f, s) and the GO magnification as µ GO .

C. Numerical integration
Given a fixed source and observer, our objective is to find the connecting GSHE trajectories of the ϵs bundle.We numerically integrate the GSHE ODEs of Eq. (2.4), or the null geodesic ODEs recovered by substituting ϵ → 0, starting at coordinate time t = 0, source position x src and initial wave packet momentum p init (k) as discussed in Section II B.
The Boyer-Lindquist coordinates contain coordinate singularities at the BH horizon and the coordinate north and south poles.Therefore, we include the following premature integration termination conditions.First, the integration is terminated if a trajectory penetrates or passes sufficiently close by the BH horizon, so that its radial component satisfies where we set ∆H = 1 + 10 −4 .Second, we terminate trajectories whose polar angle does not satisfy θ tol ≤ θ ≤ π − θ tol , where θ tol = 10 −5 .Lastly, we optionally support early termination if the absolute value of the difference between the current and initial azimuthal angle ∆ϕ = |ϕ − ϕ src | satisfies ∆ϕ > max(2π − |ϕ obs − ϕ src |, |ϕ obs −ϕ src |) since such solutions correspond to ones that complete more than one complete azimuthal loop around the BH.We refer to trajectories that do not completely loop around the BH as "direct".If no early termination condition is met, we terminate the integration when the trajectory reaches the observer's radius r obs .The integrator then outputs x f and p f , the location and momentum vectors of the trajectory at that instant.Typically, for each source-observer configuration, there exist at least two bundles that directly connect the source and the observer, with additional bundles completely looping around the BH.Rays that loop multiple times around the BH are a general signature of strong gravitational fields and images formed under such conditions are also referred to as retrolensing or glory effect [93][94][95][96][97].
We quantify whether a choice of initial direction k (and thus initial momenta) leads to a trajectory intersecting with the observer by calculating the angular distance ∆σ between the observer and the integrated trajectory cos ∆σ = cos θ f cos θ obs + sin θ f sin θ obs cos ∆ϕ f , (2.20) where θ f and ϕ f are the polar and azimuthal angles of the trajectory at r obs , and ∆ϕ f = ϕ f − ϕ obs .However, we note that in the numerical implementation we use the more accurate haversine formula for small ∆σ [98].A trajectory is considered to intersect with the observer if ∆σ → 0 and concretely we enforce that ∆σ ≤ 10 −12 .Given the nature of the GSHE, the initial directions of the GSHE trajectories at neighboring ϵ should lie sufficiently close to each other (or to the initial geodesic direction).Therefore, we typically begin by solving for the initial direction at the highest value of ϵ that connects the source and observer, then we solve for the 2 nd highest value of ϵ in a restricted region of the former initial direction and repeat this process down to the smallest ϵ and the geodesic initial direction.
We first evaluate the ODEs symbolically in Mathematica [99], expressing them explicitly in the Boyer-Lindquist coordinates.
We then export the symbolic expressions to Julia [100] and use the DifferentialEquations.jl[101] along with Optim.jlpackage [102] to integrate the ODEs and optimize the initial conditions, respectively.The Jacobian in Eq. (2.18) is calculated using automatic differentiation implemented in ForwardDiff.jl[103].

D. Quantifying the time delay
We write the observer proper time of arrival of a GSHE trajectory emitted at coordinate time t = 0 belonging to the n th bundle as τ (n) GSHE (ϵ, s).We specifically denote the proper time of arrival of the geodesic belonging to the n th bundle as τ (n) GO , as it corresponds to the GO limit of infinite frequency.We note that (2.21) We will calculate the dispersive GSHE-to-geodesic time delay as Additionally, we will also explicitly investigate the birefringent delay between the right and left polarization states ∆τ (2.23) Having fixed the background Kerr metric mass M , or equivalently its Schwarzschild radius R s , ϵ is inversely proportional to the wave packet's frequency f .Therefore, the aforementioned time delays can be expressed directly as a function of f .Dimension-full units of time can be restored by multiplying the resulting expression by R s /2c.

E. Waveform modelling
Due to the frequency-and polarization-dependent observer proper time of arrival delay with respect to the GO propagation, ∆τ , the GSHE "delays" the circular basis frequency components of the original waveform.We write the circular basis frequency-domain unlensed waveform as h0 (f, s).With the notation of Eq. (2.22), the GSHE produces a frequency-domain waveform The sum runs over the different images, i.e. bundles connecting the source and observer.The exponential encodes the frequency-and polarization-dependent time delay, and the square root encodes the magnificationinduced amplitude scaling.
We generate the unlensed linear basis waveform in PyCBC [104], which can equivalently be described in the circular basis.The right and left circularly polarized basis vectors, e R and e L , can be related to the plus and cross linearly polarized basis vectors, e + and e × , as ) discussed, e.g. in Ref. [86].
As usual, a waveform can be inverse Fourier transformed into the time domain, where we use τ to denote the observer proper time.The waveform and detector sensitivity are typically described in the linearly polarized basis.In it, the detector strain is described as where F + and F x is the antenna response function to the plus and cross polarization [87].Equivalently, the detector strain can be expressed as a function of the circularly polarized waveforms upon a suitable redefinition of the antenna response function.
Beyond visually comparing the GSHE-corrected waveforms with their geodesic counterparts, we also quantify their mismatch for a single bundle connecting the source and observer.The mismatch between two waveforms is minimized over the merger time and phase.We denote the mismatch between h GO , the GO waveform related to the unlensed waveform by including the GO magnification µ GO , and h GSHE as where t c , ϕ c are the coalescence time and phase, respectively.The mismatch depends on the noise-weighted inner product between two waveforms where S is the noise spectral density amplitude that is set by choosing a GW detector.We assume the noise to be flat across all frequencies, S(f ) = 1, as was done, e.g., in Ref. [105].
For illustration, we now ignore the minimization of the mismatch assuming that the GSHE leaves the highfrequency part of the waveform -the merger -unchanged and express the mismatch of a single circular polarization component of a waveform.Furthermore, we assume that µ(f, s) = µ GO , i.e., that the magnification of the GSHE trajectories is equal to the GO magnification, which will prove to be a sufficiently good assumption.Because the GSHE correction is a phase shift in the frequency domain, we have ⟨h GSHE , h GSHE ⟩ = ⟨h GO , h GO ⟩ and where we explicitly wrote the dependence on the circular polarization state, and we define the "mixing" angle Therefore, we have that (2.32) We will demonstrate in Section III A that the frequency dependence of γ can be isolated from the relevant scaling set by the mutual position of the source and observer, thus further simplifying this expression.

III. RESULTS
Following the prescription of Section II, we search for bundles of connecting GSHE trajectories between a fixed source and an observer on the Kerr background metric.We investigate how the GSHE-induced time delay depends on the mutual position of the source and observer.We discover that in all cases the time delay can be well approximated as a frequency-dependent power law and that the signature of the GSHE is a frequency-dependent phase shift in the inspiral part of the observed waveform.
For each configuration, we find the initial directions of a bundle of trajectories by minimizing the angular distance ∆σ of Eq. (2.20).Typically, we search the range 10 −3 ≤ ϵ ≤ 10 −1 , with 30 logarithmically spaced ϵ values.Everywhere but in Section III A 3 we resort to studying only the directly connecting bundles of trajectories to simplify the interpretation.As an example, in Fig. 1 we show two directly connecting bundles.The GSHE trajectories appear as small deviations from the geodesic trajectories with fixed boundary conditions.
In Fig. 2, we plot an example dependence of ∆σ on the initial ingoing geodesic direction.We minimize ∆σ to find the initial directions that result in a connecting trajectory between a source and an observer.The empty central region indicates the initial directions that penetrate the BH horizon, delineating the BH shadow.We also overplot in Fig. 2 the GSHE initial directions upon increasing ϵ for s = ±2.If ϵ → 0 the GSHE initial direction coincides with the initial geodesic direction, otherwise it is twisted by an angle proportional to ϵ.Now we first characterize the frequency and polariza-tion dependence of the time delay on the system configuration in Section III A and then address its impact on the observed waveform in Section III B.

A. Time delay
In Fig. 3 we plot the GSHE-to-geodesic, ∆τ (ϵ, s), and the right-to-left, ∆τ R−L (ϵ), time of arrival delays for a particular source-observer configuration.We find that, independent of the mutual positions of the source and observer, both ∆τ (ϵ, s) and ∆τ R−L (ϵ) are well described by a power law.Therefore, we introduce for the dispersive GSHE-to-geodesic and birefringent right-to-left delay, respectively.In all cases, we find α ≈ 2 and α R−L ≈ 3. We note that in the former case both α and β have what will turn out to be only a weak dependence on the circular polarization state.The difference between the right and left polarization results in the subdominant, but nonzero, ∆τ R−L (ϵ) delay.The ϵ 2 dependence of the GSHE-to-geodesic delay can be understood as follows.First, the GSHE correction to the equations of motion is proportional to ϵ and, second, to reach the same observer, the GSHE initial direction must be rotated with respect to the geodesic initial direction (see the small lines in Fig. 2).The magnitude of this rotation is proportional to ϵ, therefore, altogether these two effects yield an approximate ϵ 2 dependence.The right-to-left delay is a comparison of two perturbed solutions, which produces an ϵ 3 dependence.On the other hand, the proportionality factors, β or β R−L , are set by the relative position of the source and observer and the BH spin.β also contains information on the polarization state of the GW.As shown in the left panel of Fig. 3, in the case of two directly connecting bundles, one of the bundles' GSHE trajectories (regardless of the polarization state) arrive with a positive time delay with respect to its geodesic time of arrival, while the other bundles' GSHE trajectories arrive with a negative time delay.We verify that this holds in all configurations that we tested.
We may express ∆τ explicitly as a function of frequency in dimension-full units of as with a similar expression for the right-to-left delay ∆τ R−L .Thus, the right-to-left delay is suppressed relative to the GSHE-to-geodesic delay by an additional power of 2c/(R s f ) and generally we have |β| ≫ |β R−L | (exemplified in Fig. 3).Numerically, we find that the GSHE trajectories have a "blind spot" approximately on the opposite side of the BH that cannot be reached, regardless of the initial emission direction.In other words, given a source close to the BH, there are spacetime points on a sphere of large r that cannot be reached by GSHE trajectories, while these points can be reached by geodesics.The location and size of the blind spot depend on the position of the source, ϵ (wavelength), polarization, and the BH spin.In the Schwarzschild metric, the blind spot is a cone whose size is ∼ 0.5 degrees for r src = 5 R s and ϵ = 0.1 (upper limit considered in this work).The size decreases with higher r src and lower ϵ, approaching zero when ϵ → 0 as there is no blind spot in the geodesic case.The blind spot is exactly centered on the opposite side of the BH in the Schwarzschild metric.For a source in the equatorial plane, increasing the BH spin slightly tilts the blind spot away from the equatorial plane, and its size remains approximately unchanged.We note that the presence of the blind spot is not a numerical defect and is instead a consequence of the GSHE equations.We verify this by inspecting where the GSHE trajectories intersect the farobserver sphere upon emission in all possible directions from the source and increasing the numerical accuracy.We leave a further investigation and discussion of the blind spot for future work.
We note that each of the main GSHE trajectory bundles has opposite signs of the time delay, cf.Fig. 3.The first image to be received has β > 0 (i.e.low frequencies delayed w.r.t.geodesic), while the second image has β < 0 (low frequencies advanced w.r.t.geodesic).As geodesics correspond to extrema of the time delay, we interpret this property as the first bundle being deformed by the GSHE into longer time delays, while the second bundle gets distorted in a way that decreases the travel time.This is analogous to standard lensing theory, where images form at extrema of the time-delay function.For a point lens, the first image corresponds to the absolute minimum and the second to a saddle point of the time delay.Angular deformations around the saddle point (as found in Fig. 2) drive the time delay closer to the global minimum, explaining the lower time delay associated with β < 0. The second GSHE bundle has negative parity (µ < 0), which is consistent with a saddle-point image in the point-lens analogy.The source is otherwise at (2 Rs, π/2, 0), observer at (50 Rs, θ obs , π) and a = 0.99.When both the source and observer are in the equatorial plane the right-to-left delay vanishes due to reflection symmetry.∆τgeo is nonzero and µGO remains finite when θ obs = π/2 because of the BH spin.
We now describe the dependence of the time delay on the mutual position of the source and observer and on the spin of the BH.The BH mass enters only when we relate ϵ to frequency and restore dimension-full units of time.To demonstrate the dependence, we vary the observer's polar angle θ obs and the radial distance r src of the source from the BH.We also study the directional dependence of the GSHE, wherein we keep the source fixed but calculate the delay as a function of the emission direction.Additionally, the variation of the BH spin and observer polar angle is discussed in Appendix B. In all cases, we place the observer at r obs = 50 R s after verifying that the time delay becomes approximately independent of r obs once the observer is sufficiently far away.When we plot the power law parameters describing the time delay, we include the 1σ error bars estimated by bootstrapping.Upon varying the location of the source or observer, we associate bundles by similarity in time of arrival and initial direction.

Dependence on the observer polar angle
We begin by showing the dependence of the power law parameters, describing the time delay, on θ obs in Fig. 4. We only consider direct bundles (i.e., no complete loops around the BH) indexed by n.The source is kept at (2 R s , π/2, 0), observer at (50 R s , θ obs , π) and a = 0.99.In all cases, we find near perfect agreement with the power law parameterized as in Eq. (3.1), according to α ≈ 2 and α R−L ≈ 3. The power law proportionality of the GSHE-to-geodesic delay is typically close to an order of magnitude larger than that of the right-to-left delay, in agreement with the example configuration shown in Fig. 3.While the GSHE-to-geodesic delay is maximized when both source and observer are located in the equatorial plane, the right-to-left delay is zero in such a case, because of the reflection symmetry about the equatorial plane.We numerically verify that this condition applies more generally whenever θ obs + θ src = π.Furthermore, in the bottom panels of Fig. 4 we plot ∆τ geo defined as This is the GO time of arrival difference between the geodesics of the two direct bundles indexed by n = 1, 2.
As expected, ∆τ geo is symmetric about θ obs = π/2 as the source is in the equatorial plane.In all cases, the temporal spacing of the directly connecting bundles is several orders of magnitude larger than the GSHE-induced delay within a single bundle.In the second bottom panel we show µ GO , the magnification factor of the geodesic trajectory of each of the two bundles, which shows a weak dependence on θ obs .The magnification factor is unique for each trajectory in the bundle and therefore is also a function of ϵs.However, we find that its dependence on ϵs is negligible, and therefore we only plot the geodesic magnification factor.In fact, it will turn out that in all cases considered in this work the ϵs dependence of the magnification is negligible and we may write that Similarly, we find that in all cases the ϵs dependence of the gravitational redshift, discussed in Eq. (2.12), is negligible and well described by the gravitational redshift of the geodesic trajectory.In all cases, the image from one bundle has positive parity and negative parity for the other bundle, which also consistently holds when varying θ obs .

Dependence on the source radial distance
In Fig. 5, we plot β, β R−L , ∆τ geo and µ GO when varying r src .We do not explicitly show the power law exponent.However, we verify that α ≈ 2 and α R−L ≈ 3 remain satisfied.The source is at (r src , π/2, 0), the observer is at (50 R s , 0.4π, 3π/4) and a = 0.99.We do not place the observer directly opposite the source, instead choosing ϕ obs = 3π/4.This ensures that one of the bundles completes an azimuthal angle of 3π/4, while the other 5π/4.When the source is moved further away from the BH the former will propagate directly to the observer without entering the strong-field regime of the BH, whereas the latter is forced to effectively sling by the BH.
Figure 5 shows that in the case of direct propagation, both β and β R−L decay exponentially as the trajectories do not experience strong gradients of the gravitational field, for example approximately β ∝ 10 −0.2rsrc/Rs .On the other hand, when the trajectories are forced to sling around the BH, we find that both β and β R−L tend to a constant, non-negligible value since regardless of how distant the source is, the trajectories pass close to the BH.This suggests that it is possible to place the source far away from the BH and still obtain strong GSHE corrections, provided that the trajectories pass by the BH as expected in strong lensing.
In the bottom left panel, we plot the temporal spacing of the two bundles, ∆τ geo , which is proportional to r src .In the bottom right panel, we plot the absolute value of µ GO .Just as before, the dependence of both magnification and gravitational redshift on ϵs is negligible.We previously noted that for the bundle that is forced to sling around the BH we obtain a GSHE correction that is approximately independent of r src .However, this bundle is also exponentially demagnified, as shown in Fig. 5, with approximately µ GO ∝ 10 −0.05rsrc/Rs .Since it is the square root of the magnification that scales the signal, despite the exponential demagnification, this configuration remains an interesting avenue for detecting the GSHE, as long as r src is not too large.

Directional dependence of the GSHE
We now report on the directional dependence of the time delay from the source point of view, considering trajectories that initially point towards the BH.We emit a GSHE trajectory from the source at the maximum value of ϵ in the direction parameterized by (k 2 , k 3 ), introduced in Eq. (2.15).Then we record the angular coordinates where this trajectory intersects a far origincentered sphere of radius 50 R s , setting that location as the "observer" for the above choice of initial direction.We find the remaining GSHE and geodesic trajectories that connect to the same point and form a bundle of trajectories.Starting with the maximum value of ϵ guarantees that we never fix an observer in the blind spot of any GSHE trajectories.
We characterize each bundle belonging to an initial choice of (k 2 , k 3 ) by β of the right-polarized rays in the left panel of Fig. 6 (note that the directions in this figure correspond to the initial directions of the GSHE rays with maximum ϵ = 0.01).Throughout, we keep the source at (5 R s , π/2, 0) and do not calculate the left-polarized rays, as those behave sufficiently similarly.This time, we do not eliminate the initial directions that result in trajectories that completely loop around the BH.We still have α ≈ 2, although a small fraction of the initial directions, particularly close to the BH horizon, deviate by ∼ 1%.The left panel of Fig. 6 shows a characteristic ring of initial directions that produce |β| ∼ 1, which approximately corresponds to the trajectories that are mapped to the point opposite side of the BH (more precisely, these trajectories are mapped close to the edge of the blind spot for the maximum ϵ = 0.01).The initial directions close to the BH horizon produce |β| ∼ 10, although these are extreme configurations that completely loop around the BH and are demagnified.The initial directions of the outgoing trajectories (not shown in Fig. 6) result in |β| lower than the minimum of the ingoing trajectories and therefore are of little interest for the detection of the GSHE.
Having demonstrated how |β| depends on the direction of the emission, we now study the dependence of the corresponding magnification factor.We again verify that the deviation of the magnification as a function of ϵ from its geodesic is at most ∼ 1%, although typically smaller by up to several orders of magnitude.In Fig. 7, we show µ GO as a function of the emission direction, matching Fig. 6.Additionally, in Fig. 8 we explicitly show a scatter plot of µ GO and β corresponding to the pixels in Fig. 6 and Fig. 7.The scatter plot displays two high |β| tails -one where |β| is positively correlated with µ GO and one where the correlation is negative.The former corresponds to the aforementioned outer green ring of Fig. 6 of bundles that approximately reach the point on the other side of the BH and are magnified as they converge into a smaller region.The latter is demagnified, as it consists of bundles that pass close to the BH horizon and are sensitive to the initial direction.Therefore, it is the outer green ring of Fig. 6 that comprises a promising landscape for observing the GSHE due to its high |β| and We calculate the fraction of the source celestial halfsphere of Fig. 6 that yields |β| > β min of the GSHE-togeodesic delay for the right-polarized rays as where the integral runs over the celestial half-sphere of ingoing trajectories, H(•) is the Heaviside step function defined as H(x) = 1 if x > 0 and 0 otherwise.We note that a fraction of the half-sphere is covered by the shadow of the BH and, therefore, Υ src (0) < 1.We plot Υ src (β min ) in Fig. 9 for sources at (r src , π/2, 0), where we choose r src = 5, 7.5, 10 R s and a = 0.99.We find that for r src = 5 R s about 5% of the ingoing half-sphere yield |β| ≳ 0.5, and we verify that Υ src is approximately proportional to 1/r 2 src in the region where it is decaying.Similarly, we calculate the fraction of the far sphere of radius r = r obs where an observer would measure |β| > β min and µ > |µ min |: Here, (θ, ϕ) are coordinates on the spacetime sphere r = r obs , and (ρ, ψ) are coordinates on the celestial sphere of the source.The Jacobian relating both coordinates is the inverse of the magnification, as has been included in the second line: this can be intuitively understood as magnified/demagnified trajectories being focused/spread out and therefore less/more likely.The integral is weighted by a selection function eliminating trajectories that are either too faint to be detected or for which the GSHE is undetectable.We are considering trajectories that loop around the BH.Therefore, multiple trajectories can reach an observer, so Υ obs > 1 in general when computing probabilities (Section IV E).Fig. 10 shows the observer's cumulative GSHE probability for different magnification cuts.Two cases are considered: the left panel allowing for any number of loops around the BH, which has a maximum number of 7 in our numerical exploration.The right panel restrict the results to zero loops, although strongly deflected trajectories with α < 2π are still considered (these trajectories could be discriminated by the sign of µ, as they have negative parity).The differences are noticeable only for faint trajectories with |µ| ≪ 1: for β min ≲ 1 Υ obs is larger than unity, reflecting the existence of these additional trajectories.For β min ≳ 1 the additional loops increase the probability considerably.Note that the high β end is restricted by the resolution in our numerical exploration.
The results can be adapted to different distances between the source and the BH without an additional sampling.The GSHE probability for the source scales as Υ src ∝ r −2 src , cf.Fig. 9, as the regions contributing to the different values of β span a smaller portion of the source's sphere.Additionally, the magnification scales by the same factor µ ∝ r −2 src [59], reflecting the divergence of rays before encountering the lens.
Lastly, in Appendix C we discuss the relation between the image parity of trajectory bundles of Fig. 6 and the sign of the GSHE-to-geodesic delay.Appendix D discusses the effect of multiple loops and sign of β on the observer's probability.

Dependence on the remaining parameters
We postpone the discussion of varying the BH spin a and the azimuthal angle of the observer ϕ obs to Appendix B. However, we highlight that in the Schwarzschild metric, the right-to-left delay vanishes because of reflection symmetry.On the other hand, the GSHE-to-geodesic delay is maximized in Schwarzschild, which we attribute to the fact that lowering the BH spin pushes its horizon outwards, and therefore the trajectories pass closer to the BH horizon where the gradient of the gravitational field is larger.We verify that this behavior is not a consequence of a particular sourceobserver configuration and qualitatively holds in general.

B. Waveform comparison
We consider the IMRPhenomXP waveform [106] of a 25 and 10M ⊙ binary BH merger observed at an inclination angle of 0.9π with the spin of the primary along the zaxis a z = 0.7 and 0 along the remaining axes and zero spin of the secondary.The frequency-domain waveform is generated from 40 Hz to 1024 Hz, though the merger frequency is ∼ 225 Hz.Following Eq. (2.7), we fix the background mass to achieve some maximum value ϵ max at the lower frequency limit, since ϵ ∝ 1/f .The GSHE-to-geodesic delay shifts both polarizations in approximately the same direction with respect to the geodesic, as exemplified in Figure 3. Their difference is the right-to-left delay, which is negligible in most cases.Therefore, we will focus on the difference between the GSHE-corrected and geodesic-only waveforms.In Fig. 11, we compare the right-polarization geodesiconly and GSHE-corrected waveforms for β = 2 separately if log ϵ max = −1.5, −1.This choice of β is large enough to demonstrate the GSHE, but still reasonably likely, as we showed in Fig. 6 and Fig. 9.We follow the modeling prescription of Eq. (2.24).Even in the former, more conservative ϵ max case, the effect on the waveform is clearly visible and manifested as a frequency-dependent phase shift in the inspiral phase of the merger.This is because the merger and the ringdown are propagated by higher frequency components whose GSHE correction is suppressed as ∼ 1/f 2 .Consequently, the intrinsic parameters inferred from the inspiral part of the waveform may appear inconsistent with the merger and ringdown part of the waveform if the GSHE is not taken into account.We do not explicitly show the detector strain, which is a linear combination of the right-and left-polarization state waveforms whose phase difference due to the GSHE is negligible.In Fig. 12 we plot the mismatch of the rightpolarization waveform calculated following Eq.(2.30).We assume that in ∆τ the exponent is α = 2.We show the mismatch for several choices of ϵ max , which is equivalent to scaling the background mass M while keeping the waveform fixed.Following Eq. (2.32), this shows that we can approximate the mismatch as M ∝ β 2 for small mixing angles γ.

IV. DISCUSSION
In the derivation of the GSHE and throughout this work, several simplifying assumptions have been made to demonstrate the viability of this effect for future detection.We now first comment on the neglected higherorder contributions to the GSHE in Section IV A, the source-observer placement in Section IV B and the GW emission modelling in Section IV C.Then, in Section IV E we discuss the prospects of detecting the GSHE and, finally, in Section IV D we discuss its relation to tests of GR and beyond-GR theories.

A. Higher-order GSHE contributions
The GSHE equations describe the motion of a wave packet energy centroid and are only valid up to first order in wavelength.The relevant indicator is the WKB perturbation parameter ϵ, the ratio between the wave packet wavelength and the BH Schwarzschild radius.In the limit of ϵ → 0 the geodesic propagation of the wave packet is recovered, while ϵ ∼ 1 is the regime of wavelike phenomena, wherein the wavelength is comparable to the characteristic length scale of the system.Going further, if ϵ → ∞ we do not expect wave propagation to be significantly affected by the presence of the BH as in this limit the presence of the BH becomes negligible (see, for example, Ref. [107, Fig. 2]).
The terms of order ϵ 2 and higher were neglected in the derivation of the GSHE.In this work, we use a maximum value of ϵ = 0.1, at which point we assume that the beyond-linear terms are still negligible.Nevertheless, in Fig. 12 we showed that the effect is significant even when this maximum ϵ is relaxed.The neglected higher-order contributions are likely to induce wave-like phenomena, such as diffraction, as we depart further from the regime of GO.However, the GSHE treatment describes the motion of the energy centroid of a wave packet, which is only well defined if ϵ ≪ 1.When the wavelength reaches ϵ ∼ 1 the WKB expansion up to an arbitrary order in ϵ becomes of little interest, as the perturbation series in ϵ inevitably breaks down.Therefore, instead of extending the WKB analysis to higher orders, it is potentially more instructive to directly solve the linearized gravity perturbation propagation via, e.g., the Teukolsky equation approach [107][108][109].This approach was used to study GW emission in hierarchical triple systems in Ref. [110].An alternative but no simpler route would be a path integral approach of summing over all possible paths connecting the source and observer, whose extremum would be the classical trajectories considered in this work [111].The upside of the former treatment is its validity up to an arbitrary ϵ.Moreover, it would allow matching the GSHE results in an appropriate limit.

B. Source and observer placement
We assumed that both the observer and the source are static.The assumption of a static, far observer in the Kerr metric is a good approximation if we consider r obs → ∞, as would be the case for astrophysical observations.Throughout this work, we ensured that our conclusions are independent of the distance of the observer from the BH.Additionally, one needs to consider the gravitational and cosmological redshift.We verified that the gravitational redshift due to escaping the strongfield regime of the background BH has a negligible dependence on ϵ.It affects both the geodesic and GSHE rays equally, and we do not consider it further.The cosmological redshift due to the expansion of the universe is independent of the frequency and, thus, enters as a simple multiplicative factor.
On the other hand, the assumption of a static source may break down, particularly if the source is as close to the BH as we have considered above.This depends on the distance traveled by the source while the signal is emitted over the frequency range of a given detector.The former factor depends on the orbital period of the binary around the background BH where A is the semi-major axis of the orbit.The in-band duration of the signal depends on the GW source masses and intrinsic parameters.The typical range of LVK inband source duration are 0.1 − 100 s.The static-source assumption limits the validity of our results to shorter in-band events, including the more massive mergers expected in dynamical formation scenarios and AGNs.Our framework can be applied to longer events (e.g.lighter sources such as binary neutron star mergers), but only if they orbit a sufficiently massive BH, or are located sufficiently far.Source motion also needs to be accounted for if the GSHE signature is very sensitive on the source position.This can happen in strongly aligned systems, or for trajectories that undergo a very strong deflection, such as multiple loops around the BH.The static source assumption will be severely violated by stellar mass black holes emitting in the LISA frequency band.These sources have wavelengths several orders of magnitude larger than LVK sources.They evolve very slowly in frequency and can be observed over several years [71,72], completing multiple orbits around the massive BH [58].A treatment of a moving source would require the composition of the GSHE signal across multiple time steps and accounting for the Doppler effects; see Refs.[74,76].While the very low frequency (∼ mHz) enhances the GSHE corrections, the slow frequency evolution might make a detection challenging.Moreover, at such low frequencies the perturbative expansion in ϵ may break down, necessitating a treatment in the wave optics regime, unless the background BH is sufficiently massive as described in Eq. (2.7).The mismatch M, Eq. (2.28), of the GSHEinduced corrections as a function of |β| to a single bundle of connecting trajectories for several choices of the maximum perturbation strength ϵmax and the waveform of Fig. 11.
Another potential issue is whether the binary is tidally disrupted by the background BH.This can be described by the Hills mechanism [112][113][114] and a significant perturbation occurs when the tidal force induced by the background BH is of the same order as the binary's selfgravity.This effect has been estimated in Ref. [110] for hierarchical triple systems similar to the ones we are considering.For a binary with an orbital period of 1/f , tidal effects become important when the binary is at a radius In this paper, we always consider GWs with wavelengths smaller than ϵ max = 0.1.Thus, tidal effects can be safely ignored, as they only become significant if the binary is placed at the radius of r t ≲ 0.43M , which is below the event horizon of the background BH.Thus, binary disruption only affects our results indirectly, by precluding the formation of binaries with ϵ ≫ 1, which may later evolve to the range of frequencies probed by LVK.Addressing this effect requires detailed considerations on dynamical binary formation and migration beyond the scope of this work.

C. Source emission modeling
Our analysis relies on a simplified treatment of the GW source.Here we comment on the assumptions made: quasi-circular binaries, evolution in vaccuum, and isotropic emission.
We assume a quasi-circular binary system.However, binaries formed in dynamical environments are likely to have non-negligible eccentricity due to three-body interactions [78,115,116].Even in this case, we expect the eccentricity to be distinguishable from the GSHE via its different phase evolution.The phase of eccentric binaries evolves as ∝ f −34/9 [117, Eq. 3.13], different from the dispersive GSHE, whose phase is modified by ∼ f −1 (∆t ∝ f −2 ).Eccentricity can also be distinguished by the presence of higher modes in the signal, which are not induced by the GSHE.Other environmental effects can be distinguished for similar reasons.Reference [74] computed environmental corrections for the phase ∝ f −13/3 (accretion, acceleration) and f −16/3 (dynamical friction), distinct from that of the GSHE.We also note that these environmental effects are important for low-frequency inspirals, but very suppressed for stellar-mass coalescences.
Any other effect on the emission of GWs (such as tidal interactions with the central BH) can be included in the analysis by updating the unlensed waveform h0 (f, s) in Eq. (2.24).Even if these effects are partially degenerate, the GSHE can be unambiguously discriminated through the existence of a second image with the same underlying waveform but opposite sign of β.
We have also considered an isotropic GW emitter.However, a binary merger is an anisotropic emittersimilar to an electric dipole -and the emitted power has a directional dependence (see Ref. [59] for a treatment of strong-field lensing by Schwarzschild BHs).There are two effects in which the angular dependence of the source might play a role.First, for a given ϵs-dependent set of trajectories connecting the source and observer, the initial emission direction must be rotated away from the geodesic emission direction by an angle that is approximately proportional to ϵ.This generally corresponds to an angle of ∼ 1 degrees or lower between the low-and high-frequency components of the signal.This value is well below the sensitivity to the GW intrinsic parameters, such as the orbital inclination ι.
The angular structure of the source can cause differences among the images (bundles) formed by the background BH.The multiple images may have different relative amplitude, polarization, and merger phase, depending on which angular portion of the binary is projected onto the source for each trajectory.As an example, consider the configuration shown in Fig. 1, in which the two bundles depart in opposite directions from the source.In contrast, each GSHE trajectory encompasses an angular deviation proportional to ϵ relative to the geodesic limit for that bundle.This difference is unrelated to the GSHE corrections.However, further studies that quantify the detectability of the GSHE will need to explore this effect.

D. Relation to tests of GR
If not accounted for, the GSHE might be incorrectly interpreted as a deviation from GR.In contrast, a detection favoring beyond GR physics has to be distinguished from the GSHE.Due to its frequency dependence, the GSHE mimics three tests of GR: a modified dispersion relation, constraints of the post-Newtonian parameters, and consistency of the inspiral, merger, and ringdown phases of the signal.We will focus on the modified dispersion relation, which exactly mimics the GSHE-togeodesic time delay (i.e.β) if the right-to-left delay is negligible.The connection to the other tests is not straightforward.Hence, we will focus on the modified propagation, Eq. (4.3).
The GSHE-induced delay is degenerate with a modified dispersion relation of the form in the limit |c 0 /(hf ) 2 | ≪ 1, where h is Planck's constant.This is a particular case of a generic violation of Lorentz invariance, in which a term proportional to p n is added [118][119][120].Our case (n = 0) is equivalent to a graviton mass m 2 g = c 0 > 0 if the correction has a positive sign.However, the GSHE time delay can have either sign depending on the configuration.A modified dispersion causes a frequency-dependent time delay of a GW signal where D is an effective distance to the source that coincides with the standard luminosity distance for low redshift sources [81,Eq. 56] (see also Ref. [121]).Equating Eq. (4.4) and Eq.(3.2) yields a relation between the GW propagation and GSHE parameter The GSHE-induced delay coefficient can be probed up to a factor ∝ M D. The effective distance D is related to the source's distance (see Ref. [118,Eq. 5]), which is constrained by the amplitude of the signal.In contrast, the mass M of the background BH is unknown a-priori.Measuring M would be possible if multiple signals are received, e.g. by measuring their time delay and magnification ratio.For a single signal, it might be possible to constrain M from the orbital acceleration of the binary around the background BH, cf.Eq. (4.1).Other means of constraining M may include identifying the source's environment, e.g.via an electromagnetic counterpart, or statistically, e.g.modeling the distribution of mergers around massive BHs.
The relation in Eq. (4.5) allows us to convert LVK tests of Eq. ( 4.3) into constraints on β/M .We use the full posteriors samples from the events analyzed in the third observation run [119,120].The results are shown in Table I, where he show the 90% c.l. (confidence level) for positive and negative values of c 0 , assuming a fiducial mass of 5 × 10 4 M ⊙ .We note that the LVK analyses employ a weakly informative prior on log(c 0 ), extending many orders of magnitude below the range where the data can probe Eq.(4.3).Therefore, most of the posterior samples lie in a region that is indistinguishable from GR, leading to poor sampling of the region where data is informative.An analysis with non-logarithmic priors would lead to more efficient sampling and avoid the need to treat positive and negative values of β/c 0 separately.
The key difference between a modified dispersion relation of Eq. (4.3) and the GSHE is that the former is universal: the same coefficient c 0 represents a fundamental property of gravity and modifies the waveforms of all GW events.On the contrary, the GSHE is environmental and the correction is expected to vary between events.Therefore, to constrain β from LVK bounds on anomalous GW propagation, it is necessary to use the bounds on c 0 for individual events, rather than the combined value quoted by LVK [118][119][120].Another consequence is that GW propagation tests depend on the source distance, while the GSHE does not.Therefore, the D − c 0 correlations need to be taken into account when using Eq.(4.5) to constrain β, e.g. using the full posteriors (as in Table I).
We note that the birefringent GSHE (i.e., polarizationdependent time of arrival due to β R−L ) resembles other beyond-GR effects discussed in the literature.Scalar-tensor theories with derivative couplings to curvature [122] predict that different GW (and additional) polarization states travel at different speeds on an inhomogeneous spacetime.This birefringent effect is different from ours in three respects [82]: 1) it involves a difference in the +/× polarization, rather than R-L (right-to-left), 2) it is independent of frequency, and 3) it depends on the curvature of beyond-GR fields, which can be important over astronomical scales, rather than being confined to the vicinity of a compact object.Therefore, the time delay between polarization states associated to these theories is not bounded to any specific scale, and can range from negligible to astronomical, depending on the theory and the lensing configuration.The lack of observation of birefringence in LVK data sets stringent bounds on alternative theories [123].As deviations from GR become stronger near a compact object, detecting the GSHE imprints for mergers near a massive black hole would set extremely tight bounds on such theories.
Finally, another beyond-GR birefringence effect has been studied in Ref. [83] as emerging from higher-order corrections to GR [124,125].Like the GSHE, this form of GW birefringence involves the circular polarization states and depends on frequency, although it grows with f rather than decaying like the GSHE.Moreover, it is again assumed to be a universal property of gravity, rather than an environmental, event-dependent effect.The analysis in Ref. [83] showed that all but two GW events analyzed were compatible with GR.The outliers, GW190521 and GW191109, preferred their form of birefringence over the GR prediction.However, one cannot easily interpret this preference as due to the GSHE, as a significant β R−L is unlikely and an analogue of our, typically larger, GSHEto-geodesic delay due to β, has not been included in the analysis.Unfortunately, LIGO-Virgo did not quote any results on c 0 (Eq.(4.4)) for that event.Therefore, a more  detailed analysis would be required before reaching any conclusions.

E. Detection prospects and applications
Throughout this work we considered GW sources very close to the background BH to illustrate the consequences of the GSHE on a waveform.We have focused on the case of a background BH in the range of intermediate-mass to massive of ∼ 10 5 M ⊙ .This results in reasonable values of ϵ that make the GSHE detectable for terrestrial observatories.In case of studying the detectability of the GSHE with the longer wavelength LISA-like signals, the background BH mass would have to be correspondingly increased to achieve similar values of ϵ, such as supermassive BHs.We expect that there will be a partial degeneracy between the delay proportionality factor β and the ratio between the wavelength and the background BH mass, as both control the strength of the GSHE corrections.Nevertheless, by their definition β is independent of frequency, and therefore sufficiently high-quality data should break this degeneracy.
One of the environments to produce promising signals are AGNs, whose potential is discussed, e.g., in Ref. [116,126,127].BHs (and binaries thereof) are expected to migrate radially inward and form binaries [55,128,129].This radial migration may bring the BHs as close as ∼ 6 R s to the background BH [56].Furthermore, migration traps could promote the growth of intermediate-mass BHs around AGNs [130].In addition, a population of intermediate-mass BHs is expected in globular clusters, although no clear detection is available as of today to constrain their population [131].We consider AGNs and globular clusters to be the most likely candidates to host the hierarchical triple systems we consider, although their respective binary BH populations also remain poorly constrained [132].Although we have focused mainly on BH mergers, neutron star binaries in close proximity to an AGN would be ideal to probe the GSHE, in addition to nuclear physics [133].
We find there to exist at least two favorable sourceobserver configurations that result in a strong GSHE: aligned and close-by setups.The aligned setup occurs when the source and observer are approximately on opposite sides of the background BH.We show in Fig. 6 that in this case there exists a ring of initial directions that results in |β| ≳ 1.Because such trajectories converge to a small region opposite the source, they are also magnified, which is represented by the high |β| and high magnification cluster of points on Fig. 8. Additionally, we demonstrate that in this case it is not necessary for the source to be within a few R s of the background BH.The sufficient condition is for the trajectories to pass close to the BH.In Fig. 9, we show that the fraction of these initial directions falls approximately as 1/r 2 src .This is likely to be at least partially balanced by the fact that more mergers may occur from the outer regions of the AGN or globular cluster.
The close-by setup occurs for generic source-observer placements, but requires proximity between the source and the background BH.Even if the source, BH and the observer are not aligned, there is always a strongly deflected connecting bundle that propagates very close to the background BH and thus undergoes significant GSHE corrections.In Fig. 5, we showed that the delay proportionality factor β of such bundles tends to a constant, non-negligible value even for large separations between the source and the background BH.These trajectories exist in general, but their detectability is limited by demagnification, which is significant for sources far from the background BH and/or large deflection angles.Hence, in this setup we expect the GSHE to be detectable only for sufficiently close sources, although for most observer locations.
Our scenario predicts the reception of multiple GW signals, associated with each of the bundles connecting the source and the observer.The time delay between the signals (bundles) is proportional to the mass of the background BH, and together with the relative magnification carries information about the geometry of the system.Furthermore, each image will contain GSHE corrections of different strengths.In the aligned setup, we expect the two magnified images to have only a short time delay between them.The GSHE corrections have a sizeable β, but generally each has an opposite sign, as exemplified in Fig. 3.In the close-by setup, we expect to first detect a signal with β ≪ 1, |µ| ≈ 1, followed by a demagnified one with a strong GSHE (large |β|, |µ| ≪ 1).Unless the source is very close to the background BH, the second image will likely appear as a sub-threshold trigger due to exponential demagnification.The tools developed for the search and identification of strongly lensed GWs [134,135] can be applied to searches for GSHE imprints.A possible approach to find strongly lensed GW events is to use the posterior distribution of one image as a prior for the other image, since the two should agree if they describe the same merger [136].The short time delays between signals involved in our scenario offer two advantages.First, by lowering the chance of an unrelated event being confused as another image [137] and, secondly, by narrowing down the interval within which to search for sub-threshold triggers carrying a GSHE imprint.If the signal contains higher modes, it may be possible to distinguish type II images (saddle points in the lensing potential) from type I/III (local minima/maxima) due to the lensing-induced phase shift [138][139][140].This would provide another handle on the lensing setup, as the secondary image (negative parity, lower µ) carries this phase.
The GSHE could be used to investigate the environment of GW sources.The time delay between signals associated with different bundles can be used to constrain the background BH mass M , and β can be used to infer the alignment of the source and observer and, potentially, the background BH spin.Furthermore, a detection of a nonzero β R−L would further indicate a nonzero BH spin.In addition, the source's peculiar acceleration may be used to recover information on the mass of the background BH if the static-source approximation is broken, cf.Eq. (4.1).If the acceleration can be considered constant, it will impart a ∝ f 2 correction to the phase, which can be distinguished from the GSHE.If the deviation from the static source approximation is dramatic, as expected for LISA stellar-mass sources, much more information about the orbit can be recovered, e.g.[76].
The capacity to detect GSHE corrections in GW catalogs remains largely dependent on astrophysical factors.In this exploratory work, we demonstrate that there exist plausible configurations in which the GSHE is significant.A detectability study of the GSHE would strongly depend on the prior knowledge of the background BH population, the merger rates in their environments and their location relative to the background BH.We show that the GSHE-induced mismatch can reach M ∼ 10%.Under the mismatch and signal-to-noise ratio (SNR) criterion that two waveforms are distinguishable if the product M×SNR 2 ≳ 1 [141], we expect LVK detectors to find GSHE signatures if enough stellar-mass binaries merge in the vicinity of background BHs of intermediate mass.
Recent studies of lensed gamma-ray bursts point towards a population of objects with M ∼ 10 4 M ⊙ [142][143][144], an ideal mass range to observe the GSHE.
We now estimate the prospects of GW detectors to distinguish the GSHE in a signal.To simplify the analysis, we focus on a 30 + 30M ⊙ , non-spinning, quasi-circular binary merging at a distance of r src = 5 R s from a 10 4 M ⊙ BH.We use the IMRPhenomD waveform model [145], our framework and code for detection probabilities are based on Ref. [146].We consider two setups using the LIGO (O4 curve in Ref. [147]), Cosmic Explorer (CE; [148]) and Einstein Telescope (ET; [149]) noise curves.We assume a single interferometer for simplicity: prospects will improve when considering the LVK network, multiple arm combinations in ET or a nextgeneration network of ground detectors [150] thanks to improved SNR and sky coverage.
We quantify the observational prospects by defining the effective observable volume as Here, dVz dz (z) is the comoving volume element at the source's redshift and P det (z, |µ|, ρ th ) is the fraction of signals with SNR above a given threshold.The latter depends on the ratio between the detection threshold, ρ th , the optimal SNR at the source's redshift, √ µρ opt (z), and the effect of (de)magnification is shown explicitly.The probability of observable GSHE, dΥ obs dµ (β min , |µ|), is the derivative of Eq. (3.6) with respect to |µ|.We further enforce dΥ obs d|µ| (β min , |µ|) ≤ 1, so multiple images contribute at most as one event.We include all trajectories in our analysis (excluding trajectories with multiple loops has minimal impact on results, which is dominated by strongly deflected trajectories but with with N loop = 0, cf.Fig. 10).The minimum observable value β min is determined from the mismatch (Eq.(2.28), Fig. 12) by requiring that M(β min ) > (0.327ρ opt ) −1 , where the numerical factor relates the optimal SNR to the median SNR, given P det .This threshold, known as the Lindblom criterion [151], neglects degeneracies between parameters and thus serves as a necessary condition for observability, although it may not be sufficient.
The effective observable volume, Eq. (4.6), is shown in Table II for different detectors and background BH masses.Increasing the BH mass severely reduces V G , because only strongly deflected and demagnified trajectories lead to detectable GSHE.To facilitate the interpretation, we define an effective redshift so that V c (z G ) = V G , though it should not be interpreted as a horizon.We can obtain approximate estimates of the number of detections by multiplying V G by the expected rate R of events with this characteristics (assuming it is constant) and the observation time T obs : N GSHE = V G RT obs .The probability of detection is described by a Poisson process: in the absence of GSHE signatures, the 90% limit is given by N GSHE < ln(0.1).Table II shows 90% c.l. limits on the merger rate of objects at r src = 5 R s from the background BH of different masses, assuming no GSHE detections over an observation period of 10 years.
Figure 13 illustrates the differential effective observable volume, i.e. the integrand of Eq. (4.6) for CE with binary masses of 30 + 30 M ⊙ at a source distance of 5 R s from a 10 4 M ⊙ BH.The probability is dominated by strongly deflected but demagnified trajectories, for which GSHE distortions are substantial.Highly aligned and magnified trajectories, although less likely, still contribute significantly to detections with |µ| > 1.For ET (and similarly CE), mildly demagnified trajectories can be observed up to z ∼ 1 − 10, at least if the source merges close to the background BH.
Figure 14 shows V G for different detectors, as a function of the background BH mass and the distance to the source.The scaling of probabilities and magnificatoins with r src employed is described in Sec.III A 3. The maximum redshift of the detectable region decreases as the mass of the background BH increases, since only β ≫ 1, |µ| ≪ 1 trajectories lead to observable signals.However, our estimates are constrained by the resolution of our numerical exploration.A more precise sampling of strongly bent trajectories grazing the lightring will boost the probabilities for M ≳ 10 6 M ⊙ , although detection in those cases is likely to remain difficult even for nextgeneration ground detectors.
Although the eventual detection of GSHE depends on unknown astrophysics, the above results show how prospects will improve dramatically with the nextgeneration of GW detectors.Space detectors sensitive to lower frequencies will provide a great opportunity to probe the GSHE in a different regime.LISA, operating in the mHz window, can detect stellar-mass sources years before merger, including details of their orbit against the background BH.The lower frequencies enable our per-  turbative calculations to yield distinct predictions for binaries orbiting supermassive BHs, with the caveat that orbital effects need to be included (cf.Section IV B).The GSHE will become most dramatic for a massive background BHs ∼ 10 6 M ⊙ , such as the central BH of our galaxy.Large ϵ may even allow a clear detection of leftto-right birefringence induced by the GSHE.However, treating these cases may require a non-perturbative approach (cf.Section IV A).In the future, proposed spaceborn GW detectors will provide new opportunities to search for GSHE and wave optics-induced effects on GW propagation [152][153][154][155].

V. CONCLUSION
The gravitational spin Hall effect (GSHE) describes the propagation of a polarized wave packet of finite frequencies on a background metric in the limit of a small deviation from the geometrical optics (GO) limit.We follow the GSHE prescription as presented in Refs.[24,26].There, the GSHE is derived by inserting the Wentzel-Kramers-Brillouin (WKB) ansatz into the linearized gravity action and expanding it up to first order in wavelength.The first order contributions include the spin-orbit interaction, resulting in polarization-and frequency-dependent propagation of a wave packet.GO is recovered in the limit of infinitesimal wavelength relative to the spacetime characteristic length scale, which in our work is the Schwarzschild radius of the background metric.
The results presented in this work can be framed as a fixed spatial boundary problem.We study the GSHEinduced corrections to trajectories connecting a static source and an observer as a function of frequency and polarization.In general, for a fixed source and observer, there exist at least two connecting bundles of trajectories parameterized by ϵs, with ϵ ≡ 2λ/R s and s = ±2 for gravitational waves (GWs), each of whose infinite frequency limit (ϵ → 0) is a geodesic trajectory.There exist additional bundles that loop around the background black hole (BH).Within each bundle, we compare the time of arrival of the rays as a function of ϵs with geodesic propagation.
We find that, regardless of the mutual position of the source and observer or the BH spin, the time of arrival delay follows a power law in frequency, with an exponent of 2 or 3.The former case corresponds to the dispersive GSHE-to-geodesic and the latter to the birefringent right-to-left delay.The information about the relative source-observer position and the polarization is encoded in the power law proportionality constant.The rightto-left delay is suppressed in all but the most extreme configurations, and the time delay of trajectories within a single bundle is, thus, only weakly dependent on the polarization state.Therefore, as an approximation, it can be assumed that the GSHE time of arrival is polarization independent and only a function of frequency, i.e. that the time of arrival can be parameterized by ϵ only instead of ϵs.Consequently, there is no interference between the right-and left-polarization states, as the difference is negligible for the situations we have studied.
We study the GSHE-induced time delay dependence on the relative position of the source and observer, the direction of emission and, lastly, the BH spin.We demonstrate that the GSHE predicts birefringence effects -a different time of arrival between right-and left-polarization at a fixed frequency -only on a spinning Kerr background metric.This is expected from symmetry arguments: the left and right GW polarizations are related by a parity transformation, which would leave a Schwarzschild BH invariant, but would flip the spin of a Kerr BH.
The GSHE corrections to the gravitational waveform manifest as a frequency-dependent phase shift in the inspiral phase of a waveform, the low-frequency components, whose correction is stronger.We compare an example waveform with and without the GSHE-induced delay in Fig. 11.We also calculate the GSHE-induced waveform mismatch, which can reach ∼ 10% in plausible scenarios.Without accounting for the GSHE this may be wrongly interpreted as a violation of Lorentz invariance, anomalous GW emission or an inconsistency between inspiral-merger-ringdown.Thenceforth, any detection of such an inconsistency must eliminate the GSHE before claiming the detection of new physics.
We identify two favorable configurations for detecting the GSHE.The first case, an aligned setup, closely mimics the traditional lensing scenario.In it the source and observer are approximately on opposite sides of the background BH.In this case, the fraction of initial directions that receive a significant GSHE correction falls approximately as 1/r 2 src .The second favorable configuration, a close-by source, follows from relaxing the assumption that the source and the observer are aligned with the background BH.In this case, there exist observer-source bundles of trajectories that are strongly deflected by the background BH and hence the associated signals have a strong GSHE imprint.While these signals are demagnified, they can be observed if the signal-to-noise ratio (SNR) of the source is high, it merges sufficiently close to the background BH, or both.
These scenarios can be further probed by the existence of multiple lensed signals corresponding to the different GSHE bundles.A characteristic signature is that each of the main bundles has opposite signs of the time delay: the first received signal has positive β, with low frequency components delayed relative to the geodesic.The second signal has negative β, with low frequency components advanced relative to the geodesic, in addition to a phase shift that might be detected for GW sources emitting higher harmonics [138][139][140].If current or future GW detections reveal GSHE imprints, they may be used to constrain the fraction of events near massive and intermediate-mass BHs, providing further insight into the formation channels of compact binaries.
The GSHE is distinct from other wave propagation phenomena, such as diffraction in weak-field lensing [4][5][6].These frequency-dependent modifications of the waveform are associated to lenses at cosmological distances from the source/observer, whose mass is comparable to the GW period.Besides their conceptual differences, both effects can be distinguished in data, as the GSHE time delays converge to the geodesic/GO limit as ∼ 1/f 2 , more rapidly than the 1/f of weak-field diffraction (e.g.interpreting the phase correction in Ref. [156,Eq. 11] as a time delay).
The equivalence between the frequency dependence of the GSHE and a violation of Lorentz invariance allows us to set limits using existing LIGO-Virgo-Kagra (LVK) analyses (Table I).The 90% c.l. limits can be as stringent as |β| ≲ 10 −2 , and often differ substantially for positive/negative values of the time delay.Despite potential degeneracies with other waveform parameters, these constraints are in reasonable agreement with expectations based on the mismatch with the geodesic waveform.
We then analyse detection prospects of current and proposed GW detectors on the ground.Next-generation instruments (ET, CE) have the potential to detected GSHE signatures from events near intermediate-mass BHs (M ∼ 5 × 10 4 M ⊙ ) if the merger rate within ∼ 25R s is O(1) Gpc −3 yr −1 .These estimates are conservative, as they consider a single interferometer and are limited by the resolution of our numerical studies for trajectories grazing the background BH, which dominate the probability.The sensitivity drops sharply for larger masses and separations; however, upcoming instruments in space such as LISA [71,72], TianQin and Taiji [157] in the 2030s and proposals in the decihertz [154], milihertz [153] and microhetz [155] bands offer the best prospect for observing the GSHE.Addressing the full phenomenology of the GSHE and its detectability by next-generation detectors will require extending our formalism for non-static sources and beyond the GO expansion.
We conclude that there exists potential to unambiguously detect the GSHE.This hints at an optimistic future for studying the gravitational wave propagation in strong gravitational fields, novel tests of general theory of relativity and decoding imprints of the merger environment (e.g. the spin of the lens BH if the birefringent GSHE is observable) directly from individual waveforms.The source is at (2 Rs, π/2, 0) and the observer at (50 Rs, 0.4π, π).We again have that α ≈ 2 and αR−L ≈ 3.In the Schwarzschild metric the GSHE-to-geodesic delay is maximized, while the right-to-left delay is zero.

Figure 3 .
Figure3.The dispersive GSHE-to-geodesic delay with trajectory bundles indexed by n (left panel) and the logarithm of the absolute value of the GSHE-to-geodesic delay along with the right-to-left delay for each bundle (right panel) displaying the power law dependence of the delay.The source is at (2 Rs, π/2, 0), observer at (50 Rs, 0.4π, π) and a = 0.99.

Figure 4 .
Figure 4. Time delay parametrization upon varying the polar angle of the observer θ obs .The top row shows the power law exponent of the dispersive GSHE-to-geodesic delay α and of the birefringent right-to-left delay αR−L.The middle row shows the corresponding power law proportionality factors β and βR−L.The bottom row shows the temporal spacing of the two bundles' geodesics ∆τgeo and the geodesic magnification µGO (⊕ and ⊖ indicate positive and negative parity, respectively).The source is otherwise at (2 Rs, π/2, 0), observer at (50 Rs, θ obs , π) and a = 0.99.When both the source and observer are in the equatorial plane the right-to-left delay vanishes due to reflection symmetry.∆τgeo is nonzero and µGO remains finite when θ obs = π/2 because of the BH spin.

Figure 5 .
Figure 5.Time delay parametrization upon varying the source radial distance rsrc.Similarly to Fig.4, the top row shows β and βR−L.The bottom row shows ∆τgeo and µGO.The source is otherwise at (rsrc, π/2, 0), observer at (50 Rs, 0.4π, 0.75π) and a = 0.99.We have that α ≈ 2 and αR−L ≈ 3. The n = 1 bundle completes an azimuthal angle of 5π/4 and is deflected in the strong-field regime of the BH.Consequently, β approaches a constant value, however this bundle is exponentially demagnified.

Figure 6 .Figure 7 .
Figure 6.The dispersive GSHE-to-geodesic delay parameter β as a function of the maximum ϵ = 0.01 initial momentum parameterized by k2 and k3 (Eq.(2.15)).The source is placed at (5 Rs, π/2, 0) and the "observer" is defined as the point where the ϵmax trajectory intersects a sphere of radius 50 Rs.Each pixel represents an ϵ bundle of trajectories.

Figure 8 .Figure 9 .
Figure 8.The relation between the geodesic magnification µGO and |β| in pixels of Fig. 6.The colour represents the number of complete loops N loops around the BH.The magnified region in the top right consists of the high |β| outer green ring in Fig. 6.The demagnified region in the bottom right consists of bundles that pass very close to the BH horizon.

2 µminFigure 10 .
Figure 10.Observer's cumulative GSHE probability as a function of the minimum magnification (absolute value), including all trajectories (left) and excluding trajectories that loop around the BH (right).Only trajectories with |µ| > 10 −3 are considered.Differences are appreciable only for µ ≪ 1.

Figure 11 .
Figure 11.The GSHE-corrected and geodesic-only right polarization waveforms of a 25 and 10M⊙ merger if β = 2.We show two cases of ϵmax, the perturbation strength at the lower frequency range of the waveform, along with the corresponding GSHE-induced mismatch.The GSHE manifests as a frequency-dependent phase shift in the inspiral part of the signal.
Figure 12.The mismatch M, Eq. (2.28), of the GSHEinduced corrections as a function of |β| to a single bundle of connecting trajectories for several choices of the maximum perturbation strength ϵmax and the waveform of Fig.11.

2 Figure 13 .
Figure 13.Differential effective volume, Eq. (4.6) as a function of the magnification and redshift.The plot applies to a 30 + 30M⊙ binary at 5 Rs of a 10 4 M⊙ background BH, observed by the Einstein Telescope (see text).The solid line shows the median response distance.

Figure 14 .
Figure 14.Effective volume Eq. (4.6) as a function of the background BH mass and the separation of the source.Lines show contours of equal VG for different detectors.

Table I .
90% c.l. limits on β from LVK tests of Eq. (4.3), separately for positive and negative values of c0 while assuming background BH mass of 5 × 10 4 .We also show the median total mass Mtot and luminosity distance DL.

Table II .
Effective detection volume and equivalent redshift for different detectors and background BH masses.The results assume a 30 + 30M⊙ source at a distance of rsrc = 5 Rs from the BH, with a detection threshold of ρ thr = 8.The last column displays the 10-yr 90% c.l. limits on the merger rate for events with this characteristic, assuming no observation (in units of Gpc −3 yr −1 ).