Iron Line Variability of Discoseismic Corrugation Modes
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 21 October 2018 (MN L A TEX style file v2.2)
Iron Line Variability of Discoseismic Corrugation Modes
David Tsang , (cid:63) and Iryna Butsky (cid:63) TAPIR, California Institute of Technology, M.C. 350-17, 1200 E. California Blvd., Pasadena, CA, 91125USA Physics Department, McGill University, 3600 rue University, Montreal, QC H3A 2T8, Canada
21 October 2018
ABSTRACT
Using a fast semi-analytic raytracing code, we study the variability of relativisticallybroadened Fe-K α lines due to discoseismic oscillations concentrated in the inner-mostregions of accretion discs around black holes. The corrugation mode, or c-mode, isof particular interest as its natural frequency corresponds well to the ∼ . − α line observations will allow such discoseismic models to be confirmed or ruled outas a source of particular LFQPOs. The spectral range and frequency of the variabilityof the Fe-K α line due to corrugation modes can also potentially be used to constrainthe black hole spin if observed with sufficient temporal and spectral resolution. Key words: accretion, accretion discs – hydrodynamics – waves – black hole physics– X-rays:binaries – line:profiles
Quasi-periodic oscillations (QPOs) have been detected in the rapid variability of X-ray flux from X-ray binary systems andgalactic nuclei for decades. The advent of dedicated timing instruments such as NASA’s
Rossi X-ray Timing Explorer (RXTE)(Swank 1999) allowed detailed study of both high and low frequency QPOs in galactic X-ray binaries and AGN. For X-raybinary systems, low-frequency QPOs (LFQPOs) range in frequency from ∼ . −
15 Hz, with high amplitude and coherence(
Q > ∼ − Q ∼ − α ) lines due to gravitational red shifting and Doppler shifting in black hole accretiondiscs has been used as a probe of the space time near the inner edge of the disc, with detailed observations of X-ray binaries (cid:63) Email: [email protected] (cid:13) a r X i v : . [ a s t r o - ph . H E ] J u l D. Tsang and I. Butsky m + Κ m - Κ m - j m + j r in = r ISCO m >0, n> 0 c-mode m radius fr e qu e n c y g g Figure 1.
Propagation diagram for discoseismic g-modes and c-modes. The c-modes are trapped in the propagation region between theinner disc edge at (cid:39) r ISCO and the IVR, denoted by the curve m Ω − j / Ω ⊥ , while an evanescent region exists between the IVR and theILR ( m Ω − κ ). The higher-frequency g-modes are trapped in the propagation region between the Lindblad resonances, but are stronglydamped if they encounter the corotation resonance. providing strong constraints on their black hole spins (see e.g. Beckwith & Done 2004; Fragile et al. 2005). Fe-K α lines resultfrom the fluorescence of iron atoms within the disc after absorption of incident photons from a nearby X-ray source (usuallyassumed to be either point sources or a Compton up-scattering hot corona). As the iron line acts as a probe of the innerstructure of an accretion disc, strong perturbations or tilting of the inner structure should be reflected in variability of the ironline. If these perturbations of the inner disc structure are responsible for the LFQPOs, then QPO-phase stacked observationsof the iron line emission should produce distinctive signatures, if probed with sufficient temporal and spectral resolution, suchas with ESA’s proposed LOFT mission (Feroci et al. 2012).Observationally, correlation between the broadened iron line and QPO phase has been seen in both a galactic X-raybinary, GRS 1915 + 105 using RXTE (Miller & Homan 2005), and a Seyfert 1 galaxy, NGC 3783 using XMM Newton (Tombesi et al. 2007). Variability of line profiles due to disc perturbations have previously been modelled by Karas et al.(2001) for a toy model of spiral perturbations in an accretion disc, while Schnittman et al. (2006) explored the line variabilitydue to a precessing tilted ring. While these models do provide a strong source of line variability neither provide a compellinghydrodynamic model for disc perturbation. Here we will instead focus on the line variability due to known oscillation modesof the accretion discs.In this paper we perform a detailed calculation of the iron line signatures of discoseismic corrugation mode oscillations.The presence or absence of such QPO phase dependent signatures will provide strong evidence for the viability of corrugationmodes as a source of LFQPOs. In section 2, we outline the calculation of the corrugation eigenmodes of the disc followingSilbergleit et al. (2001), (hereafter SWO). In section 3 we describe how we utilize a fast semi-analytic raytracing code we havedeveloped to calculate the disc images and iron line variability as a function of QPO phase. This code is described in greaterdetail in Appendix A. We discuss the results of our calculations in section 4, and our conclusions in section 5.
In units of G = c = M = 1 the frequencies of free-particle orbits within a relativistic accretion disc are (e.g. Okazaki et al.1987) Ω = ( r / + a ) − , (1)Ω ⊥ = Ω(1 − a/r / + 3 a /r ) / , (2) κ = Ω(1 − /r + 8 a/r / − a /r ) / , (3) c (cid:13) , 000–000 ron Line Variability a/M = 10 -3 -2 r ISCO /M = 5.996
Figure 2.
The arbitrarily scaled vertical Lagrangian displacements ( ξ z ) of the fundamental (m=1, j=1) corrugation modes with numberof radial nodes n = 0 , , a/M .For clarity, the disc images are truncatedat r/M = 22 for a/M = 10 − , r/M = 12 for a/M = 10 − , r/M = 10 for a = 0 . r/M = 6 for a/M = 0 . r/M = 3 for a/M = 0 . H/M = 0 . /
3, with inner boundary phase parameter ϑ in = π/ where a is the dimensionless black hole spin parameter, Ω is the angular velocity, Ω ⊥ is the vertical epicyclic frequency and κ is the radial epicyclic frequency. The inner edge of the disc is located at approximately the innermost stable circular orbit, r ISCO , where κ = 0.Oscillations of frequency ω , azimuthal wave number m , and vertical wave number j have critical resonant locations inthe disc. These are the Lindblad resonances , where the co-moving perturbation frequency, ˜ ω ≡ ω − m Ω, matches the radialepicyclic frequency ( κ − ˜ ω = 0); the corotation resonance where the pattern speed matches the background fluid speed(˜ ω = 0); and the vertical resonances, where the co-moving perturbation frequency matches the vertical epicyclic frequency( j Ω ⊥ − ˜ ω = 0).In the pseudo-Newtonian case, far from the resonances, the approximate WKB dispersion relation is (Okazaki et al. 1987) c s k (cid:39) ( κ − ˜ ω )( j Ω ⊥ − ˜ ω )˜ ω (4)where c s is the sound speed, and k is the radial wave number. Corrugation modes, or c-modes, are low frequency modes ofaccretion discs with vertical structure ( j (cid:54) = 0) that have propagation region between the disc inner edge, and the inner verticalresonance (IVR) (see Figure 1). Here we study the fundamental c-modes ( m = j = 1) which have oscillation frequency equalto the Lense-Thirring precession frequency at the outer edge of their propagation regions, the inner vertical resonances. Higherfrequency inertial modes, or g-modes can also be excited in the disc in the region between the Lindblad resonances, where m Ω − κ < ω < m Ω + κ , however, they are quickly damped out by absorption at the corotation resonance unless they exist inthe small region where κ peaks. Corrugation modes are also damped by interaction with the corotation resonance but at amuch smaller rate (Tsang & Lai 2009).The corrugation modes can be modelled relativistically in the Cowling approximation (no self-gravity) by utilizing theformalism of Ipser & Lindblom (1992). Here we follow the procedure of SWO in which the relativistic perturbation equationsare separated to leading order, using a separation function Ψ. As in SWO we take background disc to be a relativistic thin accretion disc (Novikov & Thorne 1973) with scale height H (cid:28) r .We take the background space-time to be Kerr with standard metric components g µν , and utilize Boyer-Linquist coordinatesunless otherwise specified. The background disc fluid four velocity is u oν = β ( t ν + Ω φ ν ) where β = r / + ar / ( r / − r / + 2 a ) / . (5)and t ν , r ν , φ ν , and θ ν are the Boyer-Lindquist coordinate unit vectors. The sound speed in the disc is given by c s = Γ / Hβ Ω ⊥ ,where Γ > c (cid:13)000
3, with inner boundary phase parameter ϑ in = π/ where a is the dimensionless black hole spin parameter, Ω is the angular velocity, Ω ⊥ is the vertical epicyclic frequency and κ is the radial epicyclic frequency. The inner edge of the disc is located at approximately the innermost stable circular orbit, r ISCO , where κ = 0.Oscillations of frequency ω , azimuthal wave number m , and vertical wave number j have critical resonant locations inthe disc. These are the Lindblad resonances , where the co-moving perturbation frequency, ˜ ω ≡ ω − m Ω, matches the radialepicyclic frequency ( κ − ˜ ω = 0); the corotation resonance where the pattern speed matches the background fluid speed(˜ ω = 0); and the vertical resonances, where the co-moving perturbation frequency matches the vertical epicyclic frequency( j Ω ⊥ − ˜ ω = 0).In the pseudo-Newtonian case, far from the resonances, the approximate WKB dispersion relation is (Okazaki et al. 1987) c s k (cid:39) ( κ − ˜ ω )( j Ω ⊥ − ˜ ω )˜ ω (4)where c s is the sound speed, and k is the radial wave number. Corrugation modes, or c-modes, are low frequency modes ofaccretion discs with vertical structure ( j (cid:54) = 0) that have propagation region between the disc inner edge, and the inner verticalresonance (IVR) (see Figure 1). Here we study the fundamental c-modes ( m = j = 1) which have oscillation frequency equalto the Lense-Thirring precession frequency at the outer edge of their propagation regions, the inner vertical resonances. Higherfrequency inertial modes, or g-modes can also be excited in the disc in the region between the Lindblad resonances, where m Ω − κ < ω < m Ω + κ , however, they are quickly damped out by absorption at the corotation resonance unless they exist inthe small region where κ peaks. Corrugation modes are also damped by interaction with the corotation resonance but at amuch smaller rate (Tsang & Lai 2009).The corrugation modes can be modelled relativistically in the Cowling approximation (no self-gravity) by utilizing theformalism of Ipser & Lindblom (1992). Here we follow the procedure of SWO in which the relativistic perturbation equationsare separated to leading order, using a separation function Ψ. As in SWO we take background disc to be a relativistic thin accretion disc (Novikov & Thorne 1973) with scale height H (cid:28) r .We take the background space-time to be Kerr with standard metric components g µν , and utilize Boyer-Linquist coordinatesunless otherwise specified. The background disc fluid four velocity is u oν = β ( t ν + Ω φ ν ) where β = r / + ar / ( r / − r / + 2 a ) / . (5)and t ν , r ν , φ ν , and θ ν are the Boyer-Lindquist coordinate unit vectors. The sound speed in the disc is given by c s = Γ / Hβ Ω ⊥ ,where Γ > c (cid:13)000 , 000–000 D. Tsang and I. Butsky
Table 1.
Properties of the fundamental c-modes for variousblack hole spins.For these systems the disc scale height is setto
H/M = 0 .
01, adiabatic index Γ = 4 /
3, with inner boundaryphase parameter ϑ in = π/ a n ω (rad s − M − ) r ISCO /M r
IVR /M ∆ r/M − . . . . . . .
834 5 . . . .
651 10 . − .
230 5 . . . .
853 5 . . . .
615 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . For hydrostatic equilibrium the vertical density and pressure profiles are given by ρ = ρ o ( r )(1 − y ) / (Γ − , (6) p = p o ( r )(1 − y ) Γ / (Γ − , (7)where the y-coordinate has been defined such that y ≡ zH ( r ) (cid:114) Γ − . (8)For a thin disc the inner region will be dominated by radiation pressure (Γ = 4 / H (Novikov & Thorne 1973).Using the formalism of Ipser & Lindblom (1992), following SWO, Perez et al. (1997) and Nowak & Wagoner (1992) andassuming a barotropic disc and perturbation we can write the perturbation equations in term of the enthalpy perturbation δh = δp/ρ where p is the pressure and ρ is the fluid density. Here the form of the perturbation is given by δ ∝ exp[ i ( mφ − iωt )].Defining the Doppler shifted perturbation frequency to be ˜ ω ≡ ω − m Ω, we follow SWO and use the convenient variable δV ≡ δhβ ˜ ω = V r ( r ) V y ( r, y ) , (9)where V y varies only slowly with r . The coupled ordinary differential equations for the perturbation are then(1 − y ) d V y dy − y Γ − dV y dy + 2˜ ω ∗ Γ − (cid:20) − ˜ ω ∗ ˜ ω ∗ (1 − y ) (cid:21) V y = 0 , (10) d V r dr − ω − κ ddr (˜ ω − κ ) dV r dr + α (˜ ω − κ ) (cid:18) − Ψ˜ ω ∗ (cid:19) V r = 0 , (11)where ˜ ω ∗ ≡ ˜ ω/ Ω ⊥ , α ≡ β g rr /c s , and Ψ is the separation function which varies slowly in r, with Ψ (cid:39) ω ∗ in the propagationregion for the c-mode. The WKB solution given by SWO is valid to order (cid:15) ≡ − ω − m (Ω − Ω ⊥ ) m Ω ⊥ . (12)For the corrugation modes the solution to the vertical equation (10) is to lowest order V y ∝ C λj ( y ) (13)where C λj (y) is the Gegenbauer polynomial (Abramowitz & Stegun 1965), with λ = (3 − Γ) / − m = j (Γ −
1) + j (3 − Γ), which allows only particular m and j to solutions to exist for a given Γ. Only the c (cid:13) , 000–000 ron Line Variability -15 -10 -5 0 5 10 15 . . . . . . . . . n = n = n = . . . . . . . . n = n = n = . . . . . . . . n = n = n = ˆ β ˆ α ˆ α ˆ α -15 -10 -5 -15 -10 -5 Figure 3.
Constant observed frequency ratio ( g = 1 / [1 + z ]) contours for impact parameters (ˆ α, ˆ β ) for a disc around a black hole withspin a/M = 0.001 and inclination µ o = 0 . , . , . r IVR , forn = 0, 1, 2 perturbations. The variability of the corrugation modes is confined between the r ISCO and r IVR , thus the spectral variationis confined to the redshift bins which have contours within the mode propagation region. axisymmetric m = 0 = j mode and the fundamental corrugation mode m = 1 = j exist for all Γ. For the m = 1 fundamentalcorrugation modes we have V y ∝ y , with separation function (to linear order) given byΨ / ˜ ω ∗ (cid:39) − (cid:15)χ (14)where χ = 3Γ − τ , which is defined such that dτdr = ˜ ω − κ (15)with τ ( r ISCO ) = 0. The WKB solution is then V r ∝ Q − / ( τ ) cos(Φ( τ ) + Φ i ) (16)where Q ( τ ) = χ (cid:15)α ˜ ω − κ , (17)Φ i is determined by the inner boundary condition, andΦ( τ ) = (cid:90) τ Q / ( τ (cid:48) ) dτ (cid:48) . (18)For the inner boundary condition we utilize a parameterization ddr V r ( r ISCO ) cos ϑ in − V r ( r ISCO ) sin ϑ in = 0 , (19)where ϑ in generalizes our ignorance of the boundary condition at ISCO. SWO carefully analyze the singularity at ISCOas the sound speed vanishes, c s ( r ) →
0, which occurs when the torque vanishes at the inner-disc edge, as in the standardNovikov-Thorne model (Novikov & Thorne 1973). However, recent MHD simulations of accretion discs (Noble et al. 2010;Penna et al. 2010, 2012), show that for realistic discs including magnetic effects, the torque at ISCO is finite, with c s → h →
0, due to magnetic stresses connecting the ISCO material to the material in the plunging region. This provides anon-zero sound speed, resulting in well behaved non-singular behaviour at the inner boundary. In this work, for simplicity, wearbitrarily select the inner boundary ϑ in = π/ V r ( r ISCO ) = 0.With this parameterization we find tan Φ i = d ln Q/dr − tan( ϑ in ) Q / (˜ ω − κ ) . (20)Using asymptotic matching (see e.g. Tsang & Lai 2008) across the IVR, and taking n to be the number of radial nodes in thetrapping region, we can then write the WKB eigenvalue condition,Φ i + (cid:90) τ IVR Q / ( τ (cid:48) ) dτ (cid:48) = π ( n − / , (21)which can be solved for the approximate mode frequency.Using this WKB estimate for the mode frequency as an initial guess, along with the lowest order vertical solution and c (cid:13)000
0, due to magnetic stresses connecting the ISCO material to the material in the plunging region. This provides anon-zero sound speed, resulting in well behaved non-singular behaviour at the inner boundary. In this work, for simplicity, wearbitrarily select the inner boundary ϑ in = π/ V r ( r ISCO ) = 0.With this parameterization we find tan Φ i = d ln Q/dr − tan( ϑ in ) Q / (˜ ω − κ ) . (20)Using asymptotic matching (see e.g. Tsang & Lai 2008) across the IVR, and taking n to be the number of radial nodes in thetrapping region, we can then write the WKB eigenvalue condition,Φ i + (cid:90) τ IVR Q / ( τ (cid:48) ) dτ (cid:48) = π ( n − / , (21)which can be solved for the approximate mode frequency.Using this WKB estimate for the mode frequency as an initial guess, along with the lowest order vertical solution and c (cid:13)000 , 000–000 D. Tsang and I. Butsky -5 0 5-505 . . . . . . . . . n = -5 0 5 . . . . . . . . . n = -5 0 5 . . . . . . . . . n = ˆ β ˆ α ˆ α ˆ α Figure 4.
Constant observed frequency ratio ( g = 1 / [1 + z ]) contours for impact parameters (ˆ α, ˆ β ) for a disc around a black holewith spin a/M = 0.5 and inclination µ o = 0 . , . , . r IVR , for n = 2 perturbations. In this region, close to r ISCO , the gravitational redshift dominates the orbital Doppler boosting. Since thecorrugation-mode propagation region is small for high spin black holes, the redshift range with spectral variability will be correspondinglysmall. separation function we can then solve the radial differential equation (11) numerically for the eigenfrequency and eigenmodesusing standard shooting methods (e.g. Press et al. 1992). With the boundary conditions given by (19), and (cid:18) dV r dr + (cid:112) − Q (˜ ω − κ ) V r (cid:19) r = r out = 0 , (22)an evanescent decay at some point in the evanescent region, r IV R < r out < r
ILR , we find the eigenmodes shown in Figure 2and Table 1, for a disc scale height
H/M = 0 .
01 with different values of black hole spin a , and number of radial nodes n . The Lagrangian displacements (Nowak & Wagoner 1992; Perez et al. 1997) are related to the perturbations as ξ r (cid:39) − ˜ ωg rr β (˜ ω − κ ) ∂∂r δV, (23) ξ z (cid:39) − ˜ ωβ (˜ ω − N z ) (cid:18) ∂∂z δV + ρA z δV (cid:19) , (24) ξ φ (cid:39) i u t u t ˜ ω (cid:18) ∂ Ω ∂r + rν z β ∆ ∗ (cid:19) ξ r , (25)where ∆ ∗ ≡ r − r + a , ν z is the vorticity, N z is the vertical Br¨unt-Vasala frequency, and A z = βN z /∂ z p . For the barotropiccase considered here A z = N z = 0 and we have the vertical Lagrangian displacement ξ z (cid:39) − β ˜ ω ∂∂z δV. (26) To provide a comparison we also examine the line variability of a simplified tilted disc model, where we assume that the thininner disc is tilted by a small amount and precessing as a solid body with the QPO oscillation frequency ω . Here we model aprecessing slightly tilted disc simply with a vertical lagrangian perturbation given by ξ z = Ar cos( φ − ωt ) , (27)where ω is the precession frequency and we set the tangent of the tilt angle A = 0 .
01. We arbitrarily take ω to be the samefrequency as for the n = 0 c-mode for each spin a . While this toy model is not a true hydrodynamical model of the disc,similar disc tilts have been seen in simulations by Fragile et al. (2007), with the disc globally precessing rather than warpingdue to the Bardeen-Petterson effect (Bardeen & Petterson 1975). c (cid:13) , 000–000 ron Line Variability φ = φ = π φ = πφ = π . . . . . . . Observed Frequency Ratio g . . . . . . I n t e n s i t y . . . . . . . Observed Frequency Ratio g . . . . . . I n t e n s i t y φ = φ = π φ = πφ = π Figure 5.
The relative intensity (arbitrary scale) of line emission for system with black hole spin a = 0 . µ o = 0 . g ≡ ν obs /ν em = 1 / (1 + z ), for limb darkening emissivity (left) and limb brightening emissivity(right). The broadened lines are shown as a function of φ , the phase of the applied n = 2 perturbation with maximum amplitude ξ z,max (cid:39) . M located at r (cid:39) . M . To calculate the effect of such perturbations on the disc spectrum and image, we developed a new fast semi-analytic raytracingmethod (based on work from Tsang 2009), which we outline in detail in Appendix A. Solving the geodesic equations utilizingthe “Mino parameter” (Drasco & Hughes 2004), we can find the position four-vectors x ν and photon four-momenta k ν at thesurface of the the accretion disc. Here the time component of the position vector x t = ∆ t is the difference in the coordinatetime t from a particular datum value (see Appendix A4 for more detail on avoiding the divergent terms) allowing us to assessthe relative coordinate time delay observed between different photon paths. In this work we ignore secondary and higher orderimages, as these contribute relatively little to the overall flux, however this is relatively simple to include in our raytracingmethods.Along each ray the quantity I ν /ν is Lorentz invariant, where I ν is the specific intensity, and ν is the photon frequency.The total redshift of the photon is defined as 1 + z ≡ /g , where g ≡ ν obs /ν em is the observed frequency ratio (not to beconfused with the metric components g µν ), hence we have I ν obs = g I ν em along a particular ray.Covariantly the observed frequency ratio, can be given by: g = k obs µ u obs µ k em ν u em ν (28)for a distant stationary observer with four velocity u obs ν , where u em ν is the four velocity of the disc material and k em ν and k obs ν are the photon four-momenta at the disc and observer respectively.To calculate images we divide up the solid angle subtended by the black hole disc in the observer’s sky into pixels, denotedby ˆ α , the impact parameter perpendicular to the spin axis, and ˆ β , the impact parameter parallel to the spin axis (see e.g.Beckwith & Done 2004).The most important observable to be calculated for each pixel is the observed flux. We have F ν obs = (cid:90) I ν obs d Ω (29)= 1 D (cid:90) (cid:90) g I ν em d ˆ αd ˆ β (30)where D is the distance from the observer to the black hole. The emissivity of the Fe-K α fluorescence depends on a incident X-ray intensity, number density and ionization state of iron atoms. Since the disc has a vertical thermal and ionization structure,the emission may be subject to an angular dependence. We separate out the angular and radial dependence for the specificemitted intensity I ν em ( r, µ em , ν em ) = R ( r ) f ( µ em ) δ ( ν em − ν o ) (31) This code is available for download at .c (cid:13)000
The relative intensity (arbitrary scale) of line emission for system with black hole spin a = 0 . µ o = 0 . g ≡ ν obs /ν em = 1 / (1 + z ), for limb darkening emissivity (left) and limb brightening emissivity(right). The broadened lines are shown as a function of φ , the phase of the applied n = 2 perturbation with maximum amplitude ξ z,max (cid:39) . M located at r (cid:39) . M . To calculate the effect of such perturbations on the disc spectrum and image, we developed a new fast semi-analytic raytracingmethod (based on work from Tsang 2009), which we outline in detail in Appendix A. Solving the geodesic equations utilizingthe “Mino parameter” (Drasco & Hughes 2004), we can find the position four-vectors x ν and photon four-momenta k ν at thesurface of the the accretion disc. Here the time component of the position vector x t = ∆ t is the difference in the coordinatetime t from a particular datum value (see Appendix A4 for more detail on avoiding the divergent terms) allowing us to assessthe relative coordinate time delay observed between different photon paths. In this work we ignore secondary and higher orderimages, as these contribute relatively little to the overall flux, however this is relatively simple to include in our raytracingmethods.Along each ray the quantity I ν /ν is Lorentz invariant, where I ν is the specific intensity, and ν is the photon frequency.The total redshift of the photon is defined as 1 + z ≡ /g , where g ≡ ν obs /ν em is the observed frequency ratio (not to beconfused with the metric components g µν ), hence we have I ν obs = g I ν em along a particular ray.Covariantly the observed frequency ratio, can be given by: g = k obs µ u obs µ k em ν u em ν (28)for a distant stationary observer with four velocity u obs ν , where u em ν is the four velocity of the disc material and k em ν and k obs ν are the photon four-momenta at the disc and observer respectively.To calculate images we divide up the solid angle subtended by the black hole disc in the observer’s sky into pixels, denotedby ˆ α , the impact parameter perpendicular to the spin axis, and ˆ β , the impact parameter parallel to the spin axis (see e.g.Beckwith & Done 2004).The most important observable to be calculated for each pixel is the observed flux. We have F ν obs = (cid:90) I ν obs d Ω (29)= 1 D (cid:90) (cid:90) g I ν em d ˆ αd ˆ β (30)where D is the distance from the observer to the black hole. The emissivity of the Fe-K α fluorescence depends on a incident X-ray intensity, number density and ionization state of iron atoms. Since the disc has a vertical thermal and ionization structure,the emission may be subject to an angular dependence. We separate out the angular and radial dependence for the specificemitted intensity I ν em ( r, µ em , ν em ) = R ( r ) f ( µ em ) δ ( ν em − ν o ) (31) This code is available for download at .c (cid:13)000 , 000–000 D. Tsang and I. Butsky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . a = . a = . a = . a = . a = . f v a r f v a r f v a r f v a r f v a r n =0 n =1 n =2g g g µ = 0.1 Figure 6.
The fraction of the intensity ( f var ) that is emitted from the variability region for each observed frequency ratio ( g ≡ / [1 + z ])bin. With cosine of the inclination angle µ o = 0 .
1, the columns represent perturbations with different numbers of radial nodes (n=0, 1,2), while the rows represent black hole spin parameters a = 10 − , − , − , . , . f ( µ em = (1 + 2 . µ em ), with radial emissivity given by R ( r ) ∝ r − . For systems with large amplitudeperturbations we can expect fractional intensity variation on the order of f var for a given redshift bin. where µ em = cos ϑ em is the cosine of the emission angle in the local frame, and ν o is the frequency of the fluorescence linein the rest frame. We take the radial dependence of the emissivity to have the form R ( r ) ∝ r − q , assuming the standardvalue q = 3, which follows the thermal dissipation of the disc (Novikov & Thorne 1973), however steeper profiles may also beappropriate (see e.g. Svoboda et al. 2012). The angular emissivity depends on the detailed vertical thermal and ionization structure of the disc atmosphere, which isbeyond the scope of this work. For simplicity we will follow Svoboda et al. (2009), and utilize two different forms of theangular emissivity: a standard limb darkening profile f ( µ em ) = 1 + 2 . µ em , (Laor 1991), and a limb brightening profile for aplane parallel atmosphere f ( µ em ) = ln(1 + µ − ) (Haardt 1993). We do not consider an isotropic angular emissivity since forthe c-modes which are mostly incompressible perturbations, such an angular profile would result in negligible variability.The emission angle ϑ em = cos − µ em is defined as the spatial angle in the local frame between the emitted photon andthe normal to the surface of the accretion disc. The unit normal ˆ n is determined using a scalar surface function F ( x α ) andthe projected spatial derivative in the frame co-moving with the disc material n µ ≡ ∇ µ F ( x α ) + u µ u ν ∇ ν F ( x α ) , ˆ n µ ≡ n µ n ν n ν (32)where u ν is the background four velocity of the disc material. For the unperturbed case we can use the surface function F ( x ) = r cos θ = 0 which defines the surface of the equatorial plane, giving ˆ n νo = − (1 /r ) θ ν , where θ ν is the unit vector in the θ -direction at the equatorial plane. The emission angle is then given as µ em = cos ϑ em ≡ ˆ n µ k µ u ν k ν , (33)where all quantities are evaluated at the emission point of the photon at the disc surface. c (cid:13) , 000–000 ron Line Variability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . a = . a = . a = . a = . a = . f v a r f v a r f v a r f v a r f v a r n =0 n =1 n =2g g g µ = 0.5 Figure 7.
The fraction of the intensity ( f var ) that is emitted from the variability region for each observed frequency ratio ( g ≡ / [1 + z ])bin. With cosine of the inclination angle µ o = 0 .
5, the columns represent perturbations with different numbers of radial nodes (n=0, 1,2), while the rows represent black hole spin parameters a = 10 − , − , − , . , . f ( µ em = (1 + 2 . µ em ), with radial emissivity given by R ( r ) ∝ r − . For systems with large amplitudeperturbations we can expect fractional intensity variation on the order of f var for a given g -bin. To extend this to the corrugation mode we define a vertical perturbation of the surface ξ z = ξ z ( t, r, φ ). The normal vectoris now defined using the surface function F ( x α ) = r cos θ − ξ z ( t, r, φ ) = 0 (34)Since the velocity of the perturbation is small compared to the background Keplerian flow, we can continue to use theapproximation u ν (cid:39) u oν . For simplicity we will also approximate the surface displacement ξ z with the vertical Lagrangiandisplacement at the disc mid-plane. We also limit our calculations to small perturbations in the limb brightening case suchthat this approximation does not produce singular values for the angular emissivity. Figures 3 and 4 show constant redshift ( g = 1 / [1 + z ]) contours for example systems with black hole spin a = 0 .
001 and 0 . n = 0 , , a = 0 .
001 and n = 2 for a = 0 . f var emitted from the variability region r < r IVR for each spectral binfor the broadened line. All of flux at the reddest parts of the broadened line emerge from the variability region for each case,while for lower spin, the blue wing and a significant portion of the rest of the spectrum should also be variable. For high spin a = 0 .
9, we see that the variability will be confined only to the most redshifted part of the spectrum, as the variability regionis very close to the black hole where gravitational redshift strongly dominates the Doppler shift.For small amplitude perturbations spectral intensity plots were derived by binning the redshifts for each pixel, withnumber of bins dependent on the number of pixels subtended by the disc. Example spectra for various phases are shown in c (cid:13)000
9, we see that the variability will be confined only to the most redshifted part of the spectrum, as the variability regionis very close to the black hole where gravitational redshift strongly dominates the Doppler shift.For small amplitude perturbations spectral intensity plots were derived by binning the redshifts for each pixel, withnumber of bins dependent on the number of pixels subtended by the disc. Example spectra for various phases are shown in c (cid:13)000 , 000–000 D. Tsang and I. Butsky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . a = . a = . a = . a = . a = . f v a r f v a r f v a r f v a r f v a r n =0 n =1 n =2g g g µ = 0.7 Figure 8.
The fraction of the intensity ( f var ) that is emitted from the variability region for each observed frequency ratio ( g ≡ / [1 + z ])bin. With cosine of the inclination angle µ o = 0 .
7, the columns represent perturbations with different numbers of radial nodes (n=0, 1,2), while the rows represent black hole spin parameters a = 10 − , − , − , . , . f ( µ em = (1 + 2 . µ em ), with radial emissivity given by R ( r ) ∝ r − . For systems with large amplitudeperturbations we can expect fractional intensity variation on the order of f var for a given g -bin. Figure 5, showing the spectral variability for different phases. This spectral variability is highlighted in Figures 9 - 13. Inthese figures for a given black hole spin parameter a , observer inclination µ o = cos θ obs , and angular emissivity law we showthe expected unperturbed spectrum of the broadened fluorescence line. We also calculate the normalized difference betweenthe broadened line spectrum at a given oscillation phase and the phase-averaged mean spectrum, the spectral variation, andplot it as a function of oscillation phase and observed frequency ratio, g , for the simple tilted disc, and the n = 0 , Using a semi-analytic relativistic ray-tracing code for the Kerr metric we have calculated the time dependent line broadeningsignature of discoseismic corrugation modes in black hole accretion discs. We have shown that detailed spectral timing of theFe-K α line variability could demonstrate if such corrugation modes are the source of LFQPOs observed in accreting blackhole systems.The spectral variability due to corrugation modes occurs in redshift ranges determined by the propagation region of themode. In contrast the simple tilted-precessing disc model shows variability across the entire spectral range of the for thebroadened line.The redshift ranges where the variability of the corrugations modes manifests can be easily understood through theg-contours of Figures 3 and 4. In Figure 3 for a = 0 .
001 the propagation regions for the c-modes r ISCO < r < r
IVR , span therange of constant g-contours. This can be directly compared to the spectral range of the variability in Figure 9. For higherspins gravitational redshift of the line can begin to dominate the Doppler boosting due to Keplerian motion. The spectralrange of the variability is reduced for a = 0 .
5, as shown in Figure 12, and drastically reduced for the a = 0 . c (cid:13) , 000–000 ron Line Variability . . . . . . a =0.001, Limb Darkening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) µ o = . µ o = . µ o = . . . . . . . a =0.001, Limb Brightening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = µ o = . µ o = . µ o = . S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) Figure 9. Upper panels : Results assuming limb darkening angular emissivity for a black hole with spin a = 0 .
001 and inclinations µ o = cos θ obs = 0 . , .
5, and 0 .
7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc r out = 20 M .The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phaseand observed frequency ratio for the simple tilted disc, and the n = 0, n = 1, and n = 2 corrugation modes. Lower panels:
Same asabove, but for the limb brightening angular emissivity.c (cid:13)000
Same asabove, but for the limb brightening angular emissivity.c (cid:13)000 , 000–000 D. Tsang and I. Butsky . . . . . . a =0.01, Limb Darkening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) µ o = . µ o = . µ o = . . . . . . . a =0.01, Limb Brightening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) µ o = . µ o = . µ o = . Figure 10. Upper panels : Results assuming limb darkening angular emissivity for a black hole with spin a = 0 .
01 and inclinations µ o = cos θ obs = 0 . , .
5, and 0 .
7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc r out = 20 M .The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phaseand observed frequency ratio for the simple tilted disc, and the n = 0, n = 1, and n = 2 corrugation modes. Lower panels:
Same asabove, but for the limb brightening angular emissivity. c (cid:13) , 000–000 ron Line Variability . . . . . . a =0.1, Limb Darkening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) µ o = . µ o = . µ o = . . . . . . . a =0.1, Limb Brightening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) µ o = . µ o = . µ o = . Figure 11. Upper panels : Results assuming limb darkening angular emissivity for a black hole with spin a = 0 . µ o = cos θ obs = 0 . , .
5, and 0 .
7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc r out = 20 M .The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phaseand observed frequency ratio for the simple tilted disc, and the n = 0, n = 1, and n = 2 corrugation modes. Lower panels:
Same asabove, but for the limb brightening angular emissivity.c (cid:13)000
Same asabove, but for the limb brightening angular emissivity.c (cid:13)000 , 000–000 D. Tsang and I. Butsky . . . . . . I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = µ o = . µ o = . µ o = . a =0.5, Limb Darkening S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) . . . . . . a =0.5, Limb Brightening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) µ o = . µ o = . µ o = . Figure 12. Upper panels : Results assuming limb darkening angular emissivity for a black hole with spin a = 0 . µ o = cos θ obs = 0 . , .
5, and 0 .
7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc r out = 20 M .The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phaseand observed frequency ratio for the simple tilted disc, and the n = 0, n = 1, and n = 2 corrugation modes. Lower panels:
Same asabove, but for the limb brightening angular emissivity. c (cid:13) , 000–000 ron Line Variability . . . . . . a =0.9, Limb Darkening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) µ o = . µ o = . µ o = . . . . . . . a =0.9, Limb Brightening I n t e n s i t y P h a s e P h a s e P h a s e P h a s e Observed Frequency Ratio g t i l t e d n = n = n = S p e c t r a l V a r i a t i o n ( N o r m a l i z e d ) µ o = . µ o = . µ o = . Figure 13. Upper panels : Results assuming limb darkening angular emissivity for a black hole with spin a = 0 . µ o = cos θ obs = 0 . , .
5, and 0 .
7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc r out = 20 M .The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phaseand observed frequency ratio for the simple tilted disc, and the n = 0, n = 1, and n = 2 corrugation modes. Lower panels:
Same asabove, but for the limb brightening angular emissivity. The patchiness of the tilted-disc spectral variation in this case is due to mainlyto resolution and redshift binning and is not physical.c (cid:13)000
Same asabove, but for the limb brightening angular emissivity. The patchiness of the tilted-disc spectral variation in this case is due to mainlyto resolution and redshift binning and is not physical.c (cid:13)000 , 000–000 D. Tsang and I. Butsky the propagation region is much smaller and strongly dominated by gravitational redshift close to the black hole. The spectralrange for variability in the a = 0 . n = 2 mode in Figure 4.While in principle the amplitude of the corrugation modes can be large, we have limited our detailed calculations tosmall amplitude perturbations as this simplifies the ray-tracing required to capture the effects on the Fe-K α emission. Ifamplitudes of the corrugation mode remain small, it would be extremely difficult to detect variability even with upcoming X-ray observatories. However, the scaled spectral dependence of the variation as a function of phase should remain qualitativelysimilar for larger amplitude oscillations. Large amplitude corrugation modes can result in order unity variation of µ em , oreven self shadowing of the oscillation region, and are expected to lead to variations on the order of ∼ f var of the flux ineach spectral bin (see Figures 6-8). While this varies significantly for different black hole spins and observer inclination,the reddest part of the iron line, where gravitational redshift dominates, will always emerge from the region closest to theblack hole where r ISCO < r < r
IVR . For bright sources, such as GRS 1915 + 105, such large variations should be detectablegiven sufficient spectral ( ∼ .
05 keV), and temporal resolution ( ∼ . a > ∼ .
5) would have frequencies above most observed LFQPOs.The signature of the radial order ( n ) of the oscillation is evident primarily in the red and blue edges of the spectralvariation. The number of nodes in the oscillation mode is reflected in the number of nodes in the red or blue wings of thespectral variation at a particular phase, as seen in Figures 9 - 13. This is particularly prominent for the low spin cases wherethe nodes are well spaced (Figures 9 - 11). However, such fine signatures in the spectral variation would be difficult to detectunless sufficiently high spectral and timing resolution observations are stacked by oscillation phase.The presence or absence of such signatures in LFQPO-phase stacked Fe-K α observations can confirm or rule out corru-gation modes as a source of LFQPOs. However, even if variability of the inner disc structure is not the source of the broadband LFQPOs, observations with high temporal and spectral resolution X-ray instruments, such as on the proposed LOFTor ATHENA missions, would allow iron line probes of inner disc structure variability. The spectral ranges and frequenciesover which this variability is seen for corrugation modes can act as an probe of black hole spin parameter complementary toexisting spin constraints from thermal and broadened iron-line observations.Our quantitative results have relied on raytracing for small amplitude oscillations. Detailed raytracing of large amplitudevertical displacements ( ξ z > ∼ H ) was beyond the scope of this paper, however the scaled spectral variability for large amplitudesshould be qualitatively similar. Finally we note that the disc models we have used have been for a simple barotropic thin disc,and should only be taken as a demonstration of physical principle. If such signatures are indeed observed, detailed modellingof the oscillation of more realistic discs and their vertical structure should be undertaken to constrain the parameters involved. ACKNOWLEDGEMENTS
DT was supported at Caltech by the Sherman Fairchild Foundation and at McGill through funding from the CanadianInstitute for Advanced Research and the Lorne Trottier Chair in Astrophysics and Cosmology. IB conducted this researchas part of a Summer Undergraduate Research Fellowship, and was supported by NASA ATP grant no. NNX11AC37G, NSFgrant no. PHY-1151197, the David and Lucile Packard Foundation, the Alfred P. Sloan Foundation, Mrs. Albert Burford, andthe Sherman Fairchild Foundation. DT would like to acknowledge helpful discussion and useful advice from Sterl Phinney,Peter Goldreich, Christian Ott, Chris Hirata, Dong Lai, Marc Favata, Anil Zenigoulu, Chad Galley, and Andrew Cumming.
APPENDIX A: RAY-TRACING
In order to calculate various spectrum and light curve properties we must first construct a simulated image of the black holeand accretion disc in the observer’s frame. In this frame the image is broken down into individual pixels of equal solid angle,and each corresponding to a single ray emitted by the accretion disc. At the observer each pixel can be indexed by the impactparameters α ( ⊥ to the spin axis projection), and β ( (cid:107) to the spin axis projection). disc Solving for the geodesics, each of theserays can be back-traced to their source, allowing us to construct a complete image of the disc as seen by a distant observer.The contravariant components of photon momenta in a Kerr metric can be given in Boyer-Linquist coordinates, assuming c (cid:13) , 000–000 ron Line Variability G = c = M BH = 1 (e.g. Misner, Thorne & Wheeler, 1973) (cid:18) dtdλ (cid:19) = ρ − (cid:20) r + a ∆ [ E ( r + a ) − L z ] − a ( aE sin θ − L z ) (cid:21) , (A1) (cid:18) drdλ (cid:19) = ρ − [( E ( r + a ) − L z a ) − ∆(( L z − aE ) + Q )] / , (A2) (cid:18) dθdλ (cid:19) = ρ − [ Q − cos θ ( L z csc θ − E a )] / , (A3) (cid:18) dφdλ (cid:19) = ρ − (cid:20) − aE + L z csc θ + a ∆ ( E ( r + a ) − L z a ) (cid:21) , (A4)where a is the black hole spin, λ is the affine parameter, and ∆ = r − r + a . E , the photon energy, L z the angularmomentum, and Q , Carter’s constant, are constants of motion.This form allows the use of simple Runge-Kutte routines to integrate out the photon paths, and are used by many authorsas a compromise between code complexity and computational speed. Care must be taken at the turning points of the u and µ variables to ensure proper integration.Utilizing elliptical integrals, and hence greater code complexity, Cunningham and Bardeen (1973) outline a quicker methodof calculating the photon trajectories using a Hamilton-Jacobi method. Here we follow a related procedure. Rearranging andswitching variables from the affine parameter to a “Mino parameter” (see e.g. Drasco & Hughes, 2004) λ (cid:48) : dλdλ (cid:48) = ρ − E weget the following coupled first order ODEs for the coordinates as a function of mino-parameter, (cid:18) dtdλ (cid:48) (cid:19) = T ( r, θ ) ≡ (cid:20) ( r + a ) ∆ − a sin θ (cid:21) + al (cid:20) − r + a ∆ (cid:21) , (A5) (cid:18) drdλ (cid:48) (cid:19) = R ( r ) ≡ ( r + a − la ) − ( r − r + a )[( l − a ) + q ] , (A6) (cid:18) dθdλ (cid:48) (cid:19) = Θ( θ ) ≡ q − l cot θ + a cos θ, (A7) (cid:18) dφdλ (cid:48) (cid:19) = Φ( r, θ ) ≡ l csc θ + a (cid:18) r + a ∆ − (cid:19) − a l ∆ , (A8)where l = L z /E , q = Q/E . The constants of motion l and q are related to the impact parameters ( α, β ) by l = − α (cid:112) − µ o and q = β + µ o ( α − a ). By utilizing the “Mino parameter” we avoid code complications involving integrating over multipleturning points in our dependent variable (see e.g. Dexter & Agol 2009). Our code performs similarly to that of Dexter & Agol(2009), which was used as a check for computational accuracy.We can solve these ODE’s in an semi-analytic fashion using the Jacobi elliptic functions and elliptic integrals. We firstsolve for the two independant variables, r ( λ (cid:48) ) and θ ( λ (cid:48) ). A1 The R equation
We define the roots of the quartic equation R ( r ) = ( r + a − al ) − ( r − r + a )[( l − a ) + q ] = 0 , (A9) r , r , r , and r (see eg. Cadez et al., 2003), allowing us to rewrite the differential equation for dr/dλ (cid:48) as drdλ (cid:48) = (cid:112) R ( r ) (A10)= − (cid:112) ( r − r )( r − r )( r − r )( r − r ) (A11)∆ λ (cid:48) = − (cid:90) rr o dr (cid:112) ( r − r )( r − r )( r − r )( r − r ) (A12)= 2 (cid:112) ( r − r )( r − r ) F (cid:18) sin − (cid:115) ( r − r )( r − r )( r − r )( r − r ) , (cid:115) ( r − r )( r − r )( r − r )( r − r ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ∞ r (A13)= 2 (cid:112) ( r − r )( r − r ) sn − (cid:18)(cid:115) ( r − r )( r − r )( r − r )( r − r ) , (cid:115) ( r − r )( r − r )( r − r )( r − r ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ∞ r , (A14)where sn is the Jacobi elliptic function and F is the Jacobi elliptic integral of the first kind. We calculate the Jacobi elliptic functions and elliptic integrals utilizing standard recurrence relations modified to work with values onpart of the complex plane, combined with the Landen’s transforms to change the complex arguments.c (cid:13)000
We define the roots of the quartic equation R ( r ) = ( r + a − al ) − ( r − r + a )[( l − a ) + q ] = 0 , (A9) r , r , r , and r (see eg. Cadez et al., 2003), allowing us to rewrite the differential equation for dr/dλ (cid:48) as drdλ (cid:48) = (cid:112) R ( r ) (A10)= − (cid:112) ( r − r )( r − r )( r − r )( r − r ) (A11)∆ λ (cid:48) = − (cid:90) rr o dr (cid:112) ( r − r )( r − r )( r − r )( r − r ) (A12)= 2 (cid:112) ( r − r )( r − r ) F (cid:18) sin − (cid:115) ( r − r )( r − r )( r − r )( r − r ) , (cid:115) ( r − r )( r − r )( r − r )( r − r ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ∞ r (A13)= 2 (cid:112) ( r − r )( r − r ) sn − (cid:18)(cid:115) ( r − r )( r − r )( r − r )( r − r ) , (cid:115) ( r − r )( r − r )( r − r )( r − r ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ∞ r , (A14)where sn is the Jacobi elliptic function and F is the Jacobi elliptic integral of the first kind. We calculate the Jacobi elliptic functions and elliptic integrals utilizing standard recurrence relations modified to work with values onpart of the complex plane, combined with the Landen’s transforms to change the complex arguments.c (cid:13)000 , 000–000 D. Tsang and I. Butsky
With the negative value of the momentum corresponding to backtraced photon decreasing in r . Solving for r (∆ λ (cid:48) ) we get r (∆ λ (cid:48) ) = r ( r − r ) − r ( r − r ) sn ( u, κ r )( r − r ) − ( r − r ) sn ( u, κ r ) , (A15)where κ r = (cid:115) ( r − r )( r − r )( r − r )( r − r ) , (A16) u = (cid:112) r − r )( r − r )∆ λ (cid:48) − u ∞ , (A17) u ∞ = sn − (cid:18)(cid:114) r − r r − r , κ r (cid:19) , (A18)even for complex values of r n .To calculate the value of ∆ λ (cid:48) ( r ) we must carefully consider any turning points that may be encountered in r . If any root r n is a positive real value greater than the horizon radius then the photon may have a turning point in r . If no turning pointis encountered before the emission point then the value of ∆ λ (cid:48) ( r ) is given by∆ λ (cid:48) ( r ) = 2 (cid:112) ( r − r )( r − r ) (cid:20) F (cid:18) sin − ψ ∞ , κ r (cid:19) − F (cid:18) sin − ψ ( r ) , κ r (cid:19)(cid:21) , (A19)where ψ ( r ) = (cid:113) ( r − r )( r − r )( r − r )( r − r ) and ψ ∞ = ψ ( r ) | r →∞ = (cid:113) ( r − r )( r − r ) .If a photon turning point is encountered by the backtrace before reaching the emission point (ie ∆ λ (cid:48) > ∆ λ (cid:48) ( r turn )) thenthe corresponding value of ∆ λ (cid:48) is given by∆ λ (cid:48) ( r ) = 2 (cid:112) ( r − r )( r − r ) (cid:20) F (cid:18) sin − ψ ( r ) , κ r (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ∞ r turn + F (cid:18) sin − ψ ( r ) , κ r (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) r em r turn (cid:21) . (A20)As the change in sign corresponds to the change in the sign of the photon momenta at the turning point. A2 The Θ equation Examining the θ equation we perform the substitution z = cos θ such that dθdλ (cid:48) = ± (cid:114) − a z − z [ q + l − a ] + q − z (A21)= ± (cid:114) a ( z + − z )( z − z − )1 − z , (A22)where z ± = − q + l − a a ± (cid:113) ( q + l − a ) a + q a If we take χ : z = z + cos χ we have dχdθ = ± (cid:115) − z ( z + − z ) , (A23)which gives dχdλ (cid:48) = dχdθ dθdλ (cid:48) = (cid:112) a ( z − z − ) (A24)= (cid:112) a ( z + cos χ − z − ) . (A25)This gives the integral λ (cid:48) ( χ ) − λ (cid:48) ( χ o ) = (cid:90) χχ o dχ (cid:112) a ( z + cos χ − z − ) (A26)= 1 a √ z + − z − (cid:20) F (cid:18) χ, (cid:114) z + z + − z − (cid:19) − F (cid:18) χ o , (cid:114) z + z + − z − (cid:19)(cid:21) , (A27)where F( ψ, k ) is the elliptic integral of the first kind. (In Abramowitz and Stegun notation this is F = F( ψ | m ), where m = k .)Inverting this equation and solving for θ we finally obtain: θ ( λ (cid:48) ) = cos − ( √ z + cn( a √ z + − z − λ (cid:48) + u θ o , κ θ )) , (A28)where u θ o = sgn( β ) cn − (cos( θ o ) / √ z + , κ θ ) , (A29) κ θ = (cid:114) z + z + − z − , (A30)and cn( u, κ ) is the Jacobi elliptic function which has inverse cn − ( u, κ ) = F(cos − ( u ) , κ ). c (cid:13) , 000–000 ron Line Variability In order to calculate the mino-parameter corresponding to a particular value of θ we see∆ λ (cid:48) ( θ ) = 1 a √ z + − z − (cid:20) F (cid:18) χ, (cid:114) z + z + − z − (cid:19) − F (cid:18) χ o , (cid:114) z + z + − z − (cid:19)(cid:21) , (A31)where χ = cos − (cid:18) cos θ √ z + (cid:19) , (A32) χ o = sgn( β ) cos − (cid:18) cos θ o √ z + (cid:19) , (A33)thus we can solve for λ (cid:48) corresponding to intersection of the ray with simple fixed θ em disc configurations. The flat disc modelof θ em = π/ A3 The Φ Equation
The φ ( λ (cid:48) ) differential equation is more complicated than the equations for the first two spatial coordinates, however thedifferential equation can be solved by breaking the integration into two parts, integration over θ and integration over r .First rewriting the Φ( r, θ ) equation we see dφdλ (cid:48) = l − cos θ + a r + lar − r + a . (A34)Integrating the first term we see ∆ φ = (cid:90) λ em λ (cid:48) o ldλ (cid:48) − cos θ (A35)= (cid:90) λ (cid:48) em λ (cid:48) o ldλ (cid:48) − z + cn ( a √ z + − z − λ (cid:48) + u θ o , k θ ) (A36)= (cid:90) λ em λ (cid:48) o ldλ (cid:48) − z + [1 − sn ( u, k θ )] (A37)= l (1 − z + ) (cid:90) λ (cid:48) em λ (cid:48) o dλ (cid:48) z + − z + sn ( u, k θ ) (A38)= l/a (1 − z + ) √ z + − z − (cid:90) u em u o du z + − z + sn ( u, k θ ) (A39)= l/a (1 − z + ) √ z + − z − Π (cid:18) χ, z + − z + , κ θ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) χ em χ o , (A40)where u = a √ z + − z − λ (cid:48) + u o , χ = am( u, κ θ ) is the Jacobi amplitude of u , and Π( ψ, n, κ ) is the elliptic integral of the thirdkind. All these values are real, and there is no difficulty in evaluating Π( ψ, n, κ ) through the standard recurrence relations.The second term proves slightly more complicated:∆ φ = a (cid:90) λ (cid:48) em λ (cid:48) o (2 r + la ) dλ (cid:48) r − r + a (A41)= a (cid:90) r em ∞ r + la ( r − r + )( r − r − ) dλ (cid:48) dr dr (A42)= a (cid:18) r + + alr + − r − (cid:19) (cid:90) ∞ r em r − r + ) dr (cid:112) ( r − r )( r − r )( r − r )( r − r ) − a (cid:18) r − + alr + − r − (cid:19) (cid:90) ∞ r em r − r − ) dr (cid:112) ( r − r )( r − r )( r − r )( r − r ) (A43)= a (cid:20) r + alr − r + a (cid:21) ∆ λ (cid:48) + 2 a (cid:112) ( r − r )( r − r ) (cid:18) r − r r + − r − (cid:19) × (cid:20) r − + al ( r − r − )( r − r − ) Π − − r + + al ( r − r + )( r − r + ) Π + (cid:21) ∞ r em (A44)where r ± = 1 ± (cid:112) − a (A45)Π ± = Π( ψ ( r ) , n ± , κ r ) (A46) n ± = ψ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) r = r ± . (A47) c (cid:13)000
The φ ( λ (cid:48) ) differential equation is more complicated than the equations for the first two spatial coordinates, however thedifferential equation can be solved by breaking the integration into two parts, integration over θ and integration over r .First rewriting the Φ( r, θ ) equation we see dφdλ (cid:48) = l − cos θ + a r + lar − r + a . (A34)Integrating the first term we see ∆ φ = (cid:90) λ em λ (cid:48) o ldλ (cid:48) − cos θ (A35)= (cid:90) λ (cid:48) em λ (cid:48) o ldλ (cid:48) − z + cn ( a √ z + − z − λ (cid:48) + u θ o , k θ ) (A36)= (cid:90) λ em λ (cid:48) o ldλ (cid:48) − z + [1 − sn ( u, k θ )] (A37)= l (1 − z + ) (cid:90) λ (cid:48) em λ (cid:48) o dλ (cid:48) z + − z + sn ( u, k θ ) (A38)= l/a (1 − z + ) √ z + − z − (cid:90) u em u o du z + − z + sn ( u, k θ ) (A39)= l/a (1 − z + ) √ z + − z − Π (cid:18) χ, z + − z + , κ θ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) χ em χ o , (A40)where u = a √ z + − z − λ (cid:48) + u o , χ = am( u, κ θ ) is the Jacobi amplitude of u , and Π( ψ, n, κ ) is the elliptic integral of the thirdkind. All these values are real, and there is no difficulty in evaluating Π( ψ, n, κ ) through the standard recurrence relations.The second term proves slightly more complicated:∆ φ = a (cid:90) λ (cid:48) em λ (cid:48) o (2 r + la ) dλ (cid:48) r − r + a (A41)= a (cid:90) r em ∞ r + la ( r − r + )( r − r − ) dλ (cid:48) dr dr (A42)= a (cid:18) r + + alr + − r − (cid:19) (cid:90) ∞ r em r − r + ) dr (cid:112) ( r − r )( r − r )( r − r )( r − r ) − a (cid:18) r − + alr + − r − (cid:19) (cid:90) ∞ r em r − r − ) dr (cid:112) ( r − r )( r − r )( r − r )( r − r ) (A43)= a (cid:20) r + alr − r + a (cid:21) ∆ λ (cid:48) + 2 a (cid:112) ( r − r )( r − r ) (cid:18) r − r r + − r − (cid:19) × (cid:20) r − + al ( r − r − )( r − r − ) Π − − r + + al ( r − r + )( r − r + ) Π + (cid:21) ∞ r em (A44)where r ± = 1 ± (cid:112) − a (A45)Π ± = Π( ψ ( r ) , n ± , κ r ) (A46) n ± = ψ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) r = r ± . (A47) c (cid:13)000 , 000–000 D. Tsang and I. Butsky
Thus we have ∆ φ (∆ λ (cid:48) ) = ∆ φ + ∆ φ . (A48)Note that when a turning point is encountered in r, the ∆ φ must be evaluated with the appropriate sign change as for ∆ λ (cid:48) ( r ).In general the parameters of the elliptical integrals of the third kind will be complex. This limits evaluation of ∆ φ usingrecurrence relations to a particular domain of complex parameter space. Evaluations outside of this domain can be performedusing numerical quadrature. A4 The T Equation
The solution to the time coordinate is the most complicated of the four Boyer-Linquist coordinates. Simplifying the T ( r, θ )equation we have: dtdλ (cid:48) = − a sin θ + ( r + a ) + 2 alrr − r + a . (A49)The first term can be integrated in a relatively straightforward fashion as:∆ t = (cid:90) a (1 + cos θ ) dλ (cid:48) (A50)= a ∆ λ (cid:48) − (cid:90) a z + cn ( u θ , κ θ )) du θ a √ z + − z − (A51)= a ∆ λ (cid:48) − az + √ z + − z − κ θ (cid:20) E(am( u θ ) , κ θ ) − (1 − κ θ ) u θ (cid:21) u em u o (A52)= a (cid:20) a + z + (1 − κ θ ) κ θ (cid:21) ∆ λ (cid:48) − az + κ θ √ z + − z − E( χ, κ θ ) (cid:12)(cid:12)(cid:12)(cid:12) χ em χ o . (A53)For the second term,∆ t ≡ ( r + a ) +2 alrr − r + a , the analytic solution can be determined as a very long combination of ellipticalintegrals of the first, second and third kinds as well as the Jacobi elliptic functions.However, the values of observer time elapsed for the intervals between photon emission and observation at infinity arenecessarily divergent. Though one could simply evaluate ∆ t only up to a large arbitrary value of r , it is better to insteadsubtract off the same infinite constant for each ray, as our interest is limited to the elapsed observer time difference betweendifferent rays, by considering the Kerr time.Remembering that the time in Kerr-coordinates is given by a transformation: dt Kerr = dt BL + (cid:18) r + a r − r + a (cid:19) dr (A54)we can then subtract off the same constant for each ray by taking∆ t (cid:48) = ∆ t Kerr − (cid:90) r bitrary r em r + a r − r + a dr (A55)= ∆ t BL + (cid:90) r em ∞ r + a r − r + a dr − (cid:90) r bitrary r em r + a r − r + a dr (A56)= ∆ t + ∆ t − (cid:90) ∞ r bitrary r + a r − r + a dr (A57)= ∆ t + (cid:90) r bitrary r em ( r + a ) + 2 alrr − r + a dr (cid:112) R ( r )+ (cid:90) ∞ r bitrary (cid:20) ( r + a ) + 2 alrr − r + a (cid:112) R ( r ) − r + a r − r + a (cid:21) dr, (A58)where r bitrary is some r beyond the disc range.These integrations are best performed with numerical quadrature as evaluation of the elliptic integrals and Jacobi ellipticfunctions prove less efficient than the numerical integration. In addition the 1 / (cid:112) R ( r ) term makes it difficult to remove thedivergent components of the analytic solution so it may be evaluated numerically. REFERENCES
Abramowitz M., Stegun I. A., 1965, Handbook of Mathematical Functions. Dover Publications, NYBardeen J. M., Petterson J. A., 1975, Astrophys. J. Lett., 195, L65 c (cid:13) , 000–000 ron Line Variability Beckwith K., Done C., 2004, Mon. Not. Roy. Astron. Soc. , 352, 353Dexter J., Agol E., 2009, Astrophys. J., 696, 1616Drasco S., Hughes S. A., 2004, Phys. Rev. D., 69, 044015Feroci M. et al., 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8443, Societyof Photo-Optical Instrumentation Engineers (SPIE) Conference SeriesFragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007, Astrophys. J., 668, 417Fragile P. C., Miller W. A., Vandernoot E., 2005, Astrophys. J., 635, 157Haardt F., 1993, Astrophys. J., 413, 680Ipser J. R., Lindblom L., 1992, Astrophys. J., 389, 392Karas V., Martocchia A., Subr L., 2001, Pub. Astron. Soc. Jap., 53, 189Kato S., 2001, Pub. Astron. Soc. Jap., 53, 1Kato S., Fukue J., 1980, Pub. Astron. Soc. Jap., 32, 377Lai D., Tsang D., 2009, Mon. Not. Roy. Astron. Soc. , 393, 979Laor A., 1991, Astrophys. J., 376, 90Machida M., Matsumoto R., 2008, Pub. Astron. Soc. Jap., 60, 613Miller J. M., Homan J., 2005, Astrophys. J. Lett., 618, L107Noble S. C., Krolik J. H., Hawley J. F., 2010, Astrophys. J., 711, 959Novikov I. D., Thorne K. S., 1973, in Black Holes (Les Astres Occlus), Dewitt C., Dewitt B. S., eds., pp. 343–450Nowak M. A., Wagoner R. V., 1991, Astrophys. J., 378, 656Nowak M. A., Wagoner R. V., 1992, Astrophys. J., 393, 697Okazaki A. T., Kato S., Fukue J., 1987, Pub. Astron. Soc. Jap., 39, 457Penna R. F., McKinney J. C., Narayan R., Tchekhovskoy A., Shafee R., McClintock J. E., 2010, Mon. Not. Roy. Astron.Soc. , 408, 752Penna R. F., S¨a Dowski A., McKinney J. C., 2012, Mon. Not. Roy. Astron. Soc. , 420, 684Perez C. A., Silbergleit A. S., Wagoner R. V., Lehr D. E., 1997, Astrophys. J., 476, 589Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in C. Cambridge University PressSchnittman J. D., Homan J., Miller J. M., 2006, Astrophys. J., 642, 420Silbergleit A. S., Wagoner R. V., Ortega-Rodr´ıguez M., 2001, Astrophys. J., 548, 335Sobczak G. J., McClintock J. E., Remillard R. A., Cui W., Levine A. M., Morgan E. H., Orosz J. A., Bailyn C. D., 2000,Astrophys. J., 531, 537Svoboda J., Dovˇciak M., Goosmann R., Karas V., 2009, Astron. Astrophys. , 507, 1Svoboda J., Dovˇciak M., Goosmann R. W., Jethwa P., Karas V., Miniutti G., Guainazzi M., 2012, ArXiv e-printsSwank J. H., 1999, Nuclear Physics B Proceedings Supplements, 69, 12Tombesi F., de Marco B., Iwasawa K., Cappi M., Dadina M., Ponti G., Miniutti G., Palumbo G. G. C., 2007, Astron.Astrophys. , 467, 1057Tsang D., Lai D., 2008, Mon. Not. Roy. Astron. Soc. , 387, 446Tsang D., Lai D., 2009, Mon. Not. Roy. Astron. Soc. , 393, 992Tsang D. C.-W., 2009, PhD thesis, Cornell UniversityVarni`ere P., Tagger M., Rodriguez J., 2012, Astron. Astrophys. , 545, A40Wagoner R. V., 1999, Phys. Rep. , 311, 259 c (cid:13)000