Search
Close
Search
Advanced Search
Search Menu

Abstract

Using a fast semi-analytic raytracing code, we study the variability of relativistically broadened Fe–Kα lines due to discoseismic oscillations concentrated in the innermost regions of accretion discs around black holes. The corrugation mode, or c-mode, is of particular interest as its natural frequency corresponds well to the ∼0.1–15 Hz range observed for low-frequency quasi-periodic oscillations (LFQPOs) for lower spins. Comparison of the oscillation phase dependent variability and quasi-periodic oscillation-phase stacked Fe–Kα line observations will allow such discoseismic models to be confirmed or ruled out as a source of particular LFQPOs. The spectral range and frequency of the variability of the Fe–Kα line due to c-modes can also potentially be used to constrain the black hole spin if observed with sufficient temporal and spectral resolution.

accretion, accretion discs, black hole physics, hydrodynamics, line: profiles, waves, X-rays: binaries

1 INTRODUCTION

Quasi-periodic oscillations (QPOs) have been detected in the rapid variability of X-ray flux from X-ray binary systems and galactic 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-ray binary systems, low-frequency quasi-periodic oscillations (LFQPOs) range in frequency from ∼0.1 to 15 Hz, with high amplitude and coherence (Q > 10), but may vary in frequency over a short (minute) time-scale. The high-frequency quasi-periodic oscillations (HFQPOs) ranging from ∼40 to 450 Hz, have stable frequencies, but low coherence (Q ∼ 2–10). While much progress has been made in the observations of QPO phenomena, the origins of QPOs remain a mystery, with many theoretical models proposed (see Lai & Tsang 2009, and references therein for review), but thus far none have been convincingly confirmed by observation.

Perhaps the most theoretically appealing model of QPO behaviour is that of discoseismic oscillation first proposed by Kato & Fukue (1980) (see Wagoner 1999; Kato 2001, for review). In these models natural oscillation modes of thin discs in general relativistic potentials are excited at the QPO frequencies. These modes are thought to modify the flux of thermal photons from the disc which are in turn Compton up-scattered to X-ray energies by a hot corona, thereby providing variability in the X-ray flux. Discoseismic models for HFQPOs have been studied extensively (e.g. Nowak & Wagoner 1991, 1992; Perez et al. 1997; Lai & Tsang 2009).

While there are several other strong candidate models for LFQPOs [e.g. the large-scale modulation of accretion flow (Sobczak et al. 2000), torus modes (e.g. Machida & Matsumoto 2008), accretion ejection instabilities (e.g. Varnière, Tagger & Rodriguez 2012)], in this paper we will focus on the observational signatures of particular discoseismic models of the LFQPOs that may be observable using existing and upcoming instruments. Corrugation modes, or c-modes, are vertical oscillation modes trapped between the inner edge of an accretion disc and the inner vertical resonance (IVR; see e.g. Tsang & Lai 2009). These modes appear as corrugations in the disc, and typically have very low eigenfrequencies matching the Lense–Thirring precision frequency at the outer edge of their trapping region. As a source of strong perturbation of the inner disc structure with eigenfrequencies in appropriate range, they compel candidates for the source of LFQPOs.

Relativistic broadening of iron (Fe–Kα) lines due to gravitational red shifting and Doppler shifting in black hole accretion discs has been used as a probe of the space–time near the inner edge of the disc, with detailed observations of X-ray binaries providing strong constraints on their black hole spins (see e.g. Beckwith & Done 2004; Fragile, Miller & Vandernoot 2005). Fe–Kα lines result from the fluorescence of iron atoms within the disc after absorption of incident photons from a nearby X-ray source (usually assumed to be either point sources or a Compton up-scattering hot corona). As the iron line acts as a probe of the inner structure of an accretion disc, strong perturbations or tilting of the inner structure should be reflected in variability of the iron line. If these perturbations of the inner disc structure are responsible for the LFQPOs, then QPO-phase stacked observations of the iron line emission should produce distinctive signatures, if probed with sufficient temporal and spectral resolution, such as with European Space Agency'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-ray binary, 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 has previously been modelled by Karas, Martocchia & Subr (2001) for a toy model of spiral perturbations in an accretion disc by Schnittman, Homan & Miller (2006) for a precessing tilted ring, and by Ingram & Done (2012) for a precessing tilted disc. While these models do provide a strong source of line variability, none provides a compelling hydrodynamic model for disc perturbation. Here we will instead focus on the line variability due to known oscillation modes of the accretion discs.

In this paper, we perform a detailed calculation of the iron line signatures of discoseismic c-mode oscillations. The presence or absence of such QPO phase dependent signatures will provide strong evidence for the viability of c-modes as a source of LFQPOs. In Section 2, we outline the calculation of the corrugation eigenmodes of the disc following Silbergleit, Wagoner & Ortega-Rodríguez (2001, hereafter SWO). In Section 3 we describe how we utilize a fast semi-analytic raytracing code we have developed to calculate the disc images and iron line variability as a function of QPO phase. This code is described in greater detail in Appendix A. We discuss the results of our calculations in Section 4 and our conclusions in section 5.

2 CORRUGATION MODES

In units of G = c = M = 1 the frequencies of free-particle orbits within a relativistic accretion disc are (e.g. Okazaki, Kato & Fukue 1987)
\begin{equation} \Omega = (r^{3/2} + a)^{-1}\,, \end{equation}
(1)
\begin{eqnarray} \Omega _\perp = \Omega (1 - 4a/r^{3/2} + 3 a^2/r^2)^{1/2}\,, \end{eqnarray}
(2)
\begin{equation} \kappa = \Omega (1-6/r + 8a/r^{3/2} - 3a^2/r^2)^{1/2}\,, \end{equation}
(3)
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, rISCO, where κ = 0.

Oscillations of frequency ω, azimuthal wavenumber m and vertical wavenumber j have critical resonant locations in the disc. These are the Lindblad resonances, where the comoving perturbation frequency, |${\tilde{\omega }}\equiv \omega - m\Omega$|⁠, matches the radial epicyclic frequency (⁠|$\kappa ^2 - {\tilde{\omega }}^2 = 0$|⁠); the corotation resonance, where the pattern speed matches the background fluid speed (⁠|${\tilde{\omega }}= 0$|⁠); and the vertical resonances, where the comoving perturbation frequency matches the vertical epicyclic frequency (⁠|$j \Omega _\perp ^2 - {\tilde{\omega }}^2 = 0$|⁠).

In the pseudo-Newtonian case, far from the resonances, the approximate WKB dispersion relation is (Okazaki et al. 1987)
\begin{equation} c_{\rm s}^2 k^2 \simeq \frac{(\kappa ^2 - {\tilde{\omega }}^2)(j\Omega _\perp ^2 - {\tilde{\omega }}^2)}{{\tilde{\omega }}^2}, \end{equation}
(4)
where cs is the sound speed and k is the radial wavenumber. c-modes, are low-frequency modes of accretion discs with vertical structure (j ≠ 0) that have propagation region between the disc inner edge and the IVR (see Fig. 1). Here we study the fundamental c-modes (m = j = 1) which have oscillation frequency equal to the Lense–Thirring precession frequency at the outer edge of their propagation regions, the IVRs. Higher frequency 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 in the small region where κ peaks. c-modes are also damped by interaction with the corotation resonance but at a much smaller rate (Tsang & Lai 2009).
Propagation diagram for discoseismic g-modes and c-modes. The c-modes are trapped in the propagation region between the inner disc edge at ≃ rISCO and the IVR, denoted by the curve mΩ − j1/2Ω⊥, while an evanescent region exists between the IVR and the Inner Lindblad resonance given by (mΩ − κ). The higher frequency g-modes are trapped in the propagation region between the Lindblad resonances, but are strongly damped if they encounter the corotation resonance.
Figure 1.

Propagation diagram for discoseismic g-modes and c-modes. The c-modes are trapped in the propagation region between the inner disc edge at ≃ rISCO and the IVR, denoted by the curve mΩ − j1/2Ω⊥, while an evanescent region exists between the IVR and the Inner Lindblad resonance given by (mΩ − κ). The higher frequency g-modes are trapped in the propagation region between the Lindblad resonances, but are strongly damped if they encounter the corotation resonance.

Open in new tabDownload slide

The c-modes can be modelled relativistically in the Cowling approximation (no self-gravity) by utilizing the formalism of Ipser & Lindblom (1992). Here we follow the procedure of SWO in which the relativistic perturbation equations are separated to leading order, using a separation function Ψ.

2.1 Basic equations

As in SWO we take background disc to be a relativistic thin accretion disc (Novikov & Thorne 1973) with scaleheight H ≪ r. We take the background space–time to be Kerr with standard metric components gμν, and utilize Boyer–Linquist coordinates unless otherwise specified. The background disc fluid four-velocity is uoν = β(tν + Ωϕν), where
\begin{equation} \beta = \frac{r^{3/2} + a}{r^{3/4}(r^{3/2} - 3r^{1/2} + 2 a)^{1/2}} \end{equation}
(5)
and tν, rν, ϕν and θν are the Boyer–Lindquist coordinate unit vectors. The sound speed in the disc is given by cs = Γ1/2HβΩ⊥, where Γ > 1 is the adiabatic index.
For hydrostatic equilibrium the vertical density and pressure profiles are given by
\begin{eqnarray} \rho = \rho _o (r) (1-y^2)^{1/(\Gamma - 1)}, \end{eqnarray}
(6)
\begin{eqnarray} p = p_o(r) (1-y^2)^{\Gamma /(\Gamma -1)}, \end{eqnarray}
(7)
where the y-coordinate has been defined such that
\begin{equation} y \equiv \frac{z}{H(r)} \sqrt{\frac{\Gamma -1}{2\Gamma }}. \end{equation}
(8)
For a thin disc the inner region will be dominated by radiation pressure (Γ = 4/3), with nearly constant scaleheight H (Novikov & Thorne 1973).
Using the formalism of Ipser & Lindblom (1992), following SWO, Perez et al. (1997) and Nowak & Wagoner (1992) and assuming a barotropic disc and perturbation we can write the perturbation equations in terms 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 |${\tilde{\omega }}\equiv \omega - m\Omega$|⁠, we follow SWO and use the convenient variable
\begin{equation} \delta V \equiv \frac{\delta h}{\beta {\tilde{\omega }}} = V_r(r) V_y(r, y), \end{equation}
(9)
where Vy varies only slowly with r. The coupled ordinary differential equations (ODEs) for the perturbation are then
\begin{eqnarray} &&{(1-y^2) \frac{{\rm d}^2V_y}{{\rm d}y^2} - \frac{2y}{\Gamma - 1}\frac{{\rm d}V_y}{{\rm d}y} + \frac{2{\tilde{\omega }}_*^2}{\Gamma -1} \left[1 + \frac{\Psi - {\tilde{\omega }}_*^2}{{\tilde{\omega }}_*^2} (1-y)^2\right] V_y}\nonumber\\ &&{\qquad \quad = 0,} \end{eqnarray}
(10)
\begin{eqnarray} \frac{{\rm d}^2 V_r}{{\rm d}r^2} &&- \frac{1}{{\tilde{\omega }}^2 - \kappa ^2} \frac{{\rm d}}{{\rm d}r}({\tilde{\omega }}^2 - \kappa ^2) \frac{{\rm d}V_r}{{\rm d}r} + \alpha ^2 ({\tilde{\omega }}^2 - \kappa ^2) \left( 1 - \frac{\Psi }{{\tilde{\omega }}_*^2}\right)V_r\nonumber\\ && = 0, \end{eqnarray}
(11)
where |${\tilde{\omega }}_* \equiv {\tilde{\omega }}/\Omega _\perp$|⁠, |$\alpha ^2 \equiv \beta ^2 g_{rr}/c_{\rm s}^2$| and Ψ is the separation function which varies slowly in r, with |$\Psi \simeq \omega _*^2$| in the propagation region for the c-mode.

2.2 WKB solution

The WKB solution given by SWO is valid to order
\begin{equation} \epsilon \equiv -\frac{\omega - m(\Omega - \Omega _\perp )}{m \Omega _\perp }. \end{equation}
(12)
For the c-modes the solution to the vertical equation (10) is to lowest order
\begin{equation} V_y \propto C_j^\lambda (y), \end{equation}
(13)
where |$C_j^\lambda$|(y) is the Gegenbauer polynomial (Abramowitz & Stegun 1965), with λ = (3 − Γ)/2(Γ − 1). This gives rise to a selection rule 2m2 = j2(Γ − 1) + j(3 − Γ), which allows only particular m and j solutions to exist for a given Γ. Only the axisymmetric m = 0 = j mode and the fundamental c-mode m2 = 1 = j exist for all Γ. For the m = 1 fundamental c-modes we have Vy ∝ y, with separation function (to linear order) given by
\begin{equation} \Psi /{\tilde{\omega }}^* \simeq 1 - \epsilon \chi _1, \end{equation}
(14)
where χ1 = 3Γ − 1.
For the radial equation the solution is given in terms of the τ, which is defined such that
\begin{equation} \frac{{\rm d}\tau }{{\rm d}r} = {\tilde{\omega }}^2 - \kappa ^2 \end{equation}
(15)
with τ(rISCO) = 0. The WKB solution is then
\begin{equation} V_r \propto Q^{-1/4}(\tau ) \cos (\Phi (\tau ) + \Phi _i), \end{equation}
(16)
where
\begin{equation} Q(\tau ) = \frac{\chi _1 \epsilon \alpha ^2}{{\tilde{\omega }}^2 - \kappa ^2}, \end{equation}
(17)
Φi is determined by the inner boundary condition, and
\begin{equation} \Phi (\tau ) = \int _0^\tau Q^{1/2} (\tau ^{\prime }) {\rm d}\tau ^{\prime }. \end{equation}
(18)
For the inner boundary condition, we utilize a parametrization
\begin{equation} \frac{{\rm d}}{{\rm d}r}V_r(r_{\rm ISCO}) \cos \vartheta _{\rm in} - V_r(r_{\rm ISCO}) \sin \vartheta _{\rm in} = 0, \end{equation}
(19)
where ϑin generalizes our ignorance of the boundary condition at innermost stable circular orbit (ISCO). SWO carefully analyse the singularity at ISCO as the sound speed vanishes, cs(r) → 0, which occurs when the torque vanishes at the inner disc edge, as in the standard Novikov–Thorne model (Novikov & Thorne 1973). However, recent magnetohydrodynamic simulations of accretion discs (Noble, Krolik & Hawley 2010; Penna et al. 2010; Penna, Sä Dowski & McKinney 2012) show that for realistic discs including magnetic effects, the torque at ISCO is finite, with cs → 0 as h → 0, due to magnetic stresses connecting the ISCO material to the material in the plunging region. This provides a non-zero sound speed, resulting in well-behaved non-singular behaviour at the inner boundary. In this work, for simplicity, we arbitrarily select the inner boundary ϑin = π/2 such that Vr(rISCO) = 0.
With this parametrization we find
\begin{equation} \tan \Phi _i= \frac{\mathrm{d} \ln Q/\mathrm{d}r - \tan (\vartheta _{\rm in})}{Q^{1/2}({\tilde{\omega }}^2 - \kappa ^2)}. \end{equation}
(20)
Using asymptotic matching (see e.g. Tsang & Lai 2008) across the IVR, and taking n to be the number of radial nodes in the trapping region, we can then write the WKB eigenvalue condition,
\begin{equation} \Phi _i + \int _0^{\tau _{\rm IVR}} Q^{1/2} (\tau ^{\prime }) d\tau ^{\prime } = \pi (n - 1/4), \end{equation}
(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 separation function we can then solve the radial differential equation (11) numerically for the eigenfrequency and eigenmodes using standard shooting methods (e.g. Press et al. 1992). With the boundary conditions given by (19), and
\begin{equation} \bigg (\frac{{\rm d}V_r}{{\rm d}r} + \sqrt{-Q}({\tilde{\omega }}^2 - \kappa ^2) V_r \bigg )_{r = r_{\rm out}}=0, \end{equation}
(22)
an evanescent decay at some point in the evanescent region, rIVR < rout < rILR, we find the eigenmodes shown in Fig. 2 and Table 1, for a disc scaleheight H/M = 0.01 with different values of black hole spin a, and number of radial nodes n.
The arbitrarily scaled vertical Lagrangian displacements (ξz) of the fundamental (m = 1, j = 1) c-modes with number of radial nodes n = 0, 1, 2 calculated following SWO, for various black hole spin parameters a/M. For clarity, the disc images are truncated at r/M = 22 for a/M = 10−3, r/M = 12 for a/M = 10−2, r/M = 10 for a = 0.1, r/M = 6 for a/M = 0.5 and r/M = 3 for a/M = 0.9, as the capture region for modes is quite small for higher spin systems. For these systems, the disc scaleheight is set to H/M = 0.01, adiabatic index Γ = 4/3, with inner boundary phase parameter ϑin = π/2.
Figure 2.

The arbitrarily scaled vertical Lagrangian displacements (ξz) of the fundamental (m = 1, j = 1) c-modes with number of radial nodes n = 0, 1, 2 calculated following SWO, for various black hole spin parameters a/M. For clarity, the disc images are truncated at r/M = 22 for a/M = 10−3, r/M = 12 for a/M = 10−2, r/M = 10 for a = 0.1, r/M = 6 for a/M = 0.5 and r/M = 3 for a/M = 0.9, as the capture region for modes is quite small for higher spin systems. For these systems, the disc scaleheight is set to H/M = 0.01, adiabatic index Γ = 4/3, with inner boundary phase parameter ϑin = π/2.

Open in new tabDownload slide
Table 1.
Open in new tab

Properties of the fundamental c-modes for various black hole spins. For these systems, the disc scaleheight is set to H/M = 0.01, adiabatic index Γ = 4/3, with inner boundary phase parameter ϑin = π/2.

an|$\omega ({\rm rad}\ {\rm s}^{-1}\, M_{10}^{-1})$|rISCO/MrIVR/MΔr/M
00.06655.99678.48802.4913
10−310.02455.996711.8345.8380
20.008815.996716.65110.654
01.2305.96736.90860.9413
10−210.8535.96737.80511.8378
20.6155.96738.70742.7401
017.85.66936.05480.3855
0.1115.25.66936.37580.7065
213.45.66936.66240.9931
01954.23304.39870.1657
0.511794.23304.53010.2971
21674.23304.64280.4098
014602.32092.38780.0669
0.9113802.32092.43940.1185
213202.32092.48280.1619
an|$\omega ({\rm rad}\ {\rm s}^{-1}\, M_{10}^{-1})$|rISCO/MrIVR/MΔr/M
00.06655.99678.48802.4913
10−310.02455.996711.8345.8380
20.008815.996716.65110.654
01.2305.96736.90860.9413
10−210.8535.96737.80511.8378
20.6155.96738.70742.7401
017.85.66936.05480.3855
0.1115.25.66936.37580.7065
213.45.66936.66240.9931
01954.23304.39870.1657
0.511794.23304.53010.2971
21674.23304.64280.4098
014602.32092.38780.0669
0.9113802.32092.43940.1185
213202.32092.48280.1619
Table 1.
Open in new tab

Properties of the fundamental c-modes for various black hole spins. For these systems, the disc scaleheight is set to H/M = 0.01, adiabatic index Γ = 4/3, with inner boundary phase parameter ϑin = π/2.

an|$\omega ({\rm rad}\ {\rm s}^{-1}\, M_{10}^{-1})$|rISCO/MrIVR/MΔr/M
00.06655.99678.48802.4913
10−310.02455.996711.8345.8380
20.008815.996716.65110.654
01.2305.96736.90860.9413
10−210.8535.96737.80511.8378
20.6155.96738.70742.7401
017.85.66936.05480.3855
0.1115.25.66936.37580.7065
213.45.66936.66240.9931
01954.23304.39870.1657
0.511794.23304.53010.2971
21674.23304.64280.4098
014602.32092.38780.0669
0.9113802.32092.43940.1185
213202.32092.48280.1619
an|$\omega ({\rm rad}\ {\rm s}^{-1}\, M_{10}^{-1})$|rISCO/MrIVR/MΔr/M
00.06655.99678.48802.4913
10−310.02455.996711.8345.8380
20.008815.996716.65110.654
01.2305.96736.90860.9413
10−210.8535.96737.80511.8378
20.6155.96738.70742.7401
017.85.66936.05480.3855
0.1115.25.66936.37580.7065
213.45.66936.66240.9931
01954.23304.39870.1657
0.511794.23304.53010.2971
21674.23304.64280.4098
014602.32092.38780.0669
0.9113802.32092.43940.1185
213202.32092.48280.1619

2.3 Lagrangian displacements

The Lagrangian displacements (Nowak & Wagoner 1992; Perez et al. 1997) are related to the perturbations as
\begin{eqnarray} \xi ^r \simeq -\frac{{\tilde{\omega }}g^{rr}}{\beta ({\tilde{\omega }}^2 - \kappa ^2)} \frac{\mathrm{\partial} }{\mathrm{\partial} r}\delta V, \end{eqnarray}
(23)
\begin{eqnarray} \xi ^z \simeq -\frac{{\tilde{\omega }}}{\beta ({\tilde{\omega }}^2 - N_z^2)} \left( \frac{\mathrm{\partial} }{\mathrm{\partial} z} \delta V + \rho A_z \delta V\right), \end{eqnarray}
(24)
\begin{eqnarray} \xi ^\phi \simeq i\frac{u^t u_t}{{\tilde{\omega }}} \left( \frac{\mathrm{\partial} \Omega }{\mathrm{\partial} r} + \frac{r \nu ^z}{\beta ^2 \Delta ^*} \right)\xi ^r, \end{eqnarray}
(25)
where Δ* ≡ r2 − 2r + a2, νz is the vorticity, Nz is the vertical Brünt–Vasala frequency and |$A_z = \beta N_z^2/\mathrm{\partial} _z p$|⁠. For the barotropic case considered here Az = Nz = 0 and we have the vertical Lagrangian displacement
\begin{equation} \xi ^z \simeq -\frac{1}{\beta {\tilde{\omega }}}\frac{\mathrm{\partial} }{\mathrm{\partial} z}\delta V. \end{equation}
(26)

2.4 Tilted discs

To provide a comparison, we also examine the line variability of a simplified tilted disc model, where we assume that the thin inner disc is tilted by a small amount and precessing as a solid body with the QPO oscillation frequency ω. Here, we model a precessing slightly tilted disc simply with a vertical Lagrangian perturbation given by
\begin{equation} \xi ^z = A r \cos (\phi - \omega t), \end{equation}
(27)
where ω is the precession frequency and we set the tangent of the tilt angle A = 0.01. We arbitrarily take ω to be the same frequency 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 warping due to the Bardeen–Petterson effect (Bardeen & Petterson 1975).

3 CALCULATING OBSERVABLES

To calculate the effect of such perturbations on the disc spectrum and image, we developed a new fast semi-analytic raytracing method1 (based on work from Tsang 2009), which we outline in detail in Appendix A. Solving the geodesic equations utilizing the ‘Mino parameter’ (Drasco & Hughes 2004), we can find the position four-vectors xν and photon four-momenta kν at the surface of the accretion disc. Here the time component of the position vector xt = Δt is the difference in the coordinate time t from a particular datum value (see Appendix A4 for more detail on avoiding the divergent terms) allowing us to assess the relative coordinate time delay observed between different photon paths. In this work, we ignore secondary and higher order images, as these contribute relatively little to the overall flux, however, this is relatively simple to include in our raytracing methods.

Along each ray the quantity Iν/ν3 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 ≡ 1/g, where g ≡ νobs/νem is the observed frequency ratio (not to be confused with the metric comp-onents gμν), hence we have |$I_{\nu _{\rm obs}} = g^3 I_{\nu _{\rm em}}$| along a particular ray.

Covariantly the observed frequency ratio, can be given by
\begin{equation} g = \frac{k^{\rm obs}{}_\mu u_{\rm obs}{}^\mu }{k^{\rm em}{}_\nu u_{\rm em}{}^\nu }, \end{equation}
(28)
for a distant stationary observer with four-velocity uobsν, where uemν is the four-velocity of the disc material and kemν and kobsν 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, denoted by |$\hat{\alpha }$|⁠, the impact parameter perpendicular to the spin axis, and |$\hat{\beta }$|⁠, 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
\begin{eqnarray} F_{\nu _{\rm obs}} = \int I_{\nu _{\rm obs}} \mathrm{d}\Omega \end{eqnarray}
(29)
\begin{eqnarray} \hphantom{F_{\nu _{\rm obs}}}= \frac{1}{D^2}\int \!\!\!\int g^3 I_{\nu _{\rm em}}\mathrm{d}\hat{\alpha } \mathrm{d}\hat{\beta }, \end{eqnarray}
(30)
where D is the distance from the observer to the black hole. The emissivity of the Fe–Kα fluorescence depends on an 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 specific emitted intensity
\begin{equation} I_{\nu _{\rm em}}(r, \mu _{\rm em}, \nu _{\rm em}) = {\cal R}(r) f(\mu _{\rm em}) \delta (\nu _{\rm em}- \nu _o), \end{equation}
(31)
where μem = cos ϑem is the cosine of the emission angle in the local frame and νo is the frequency of the fluorescence line in the rest frame. We take the radial dependence of the emissivity to have the form |${\cal R}(r) \propto r^{-q}$|⁠, assuming the standard value q = 3, which follows the thermal dissipation of the disc (Novikov & Thorne 1973), however, steeper profiles may also be appropriate (see e.g. Svoboda et al. 2012).

3.1 Angular emissivity

The angular emissivity depends on the detailed vertical thermal and ionization structure of the disc atmosphere, which is beyond the scope of this work. For simplicity, we will follow Svoboda et al. (2009), and utilize two different forms of the angular emissivity: a standard limb darkening profile f(μem) = 1 + 2.06 μem, (Laor 1991) and a limb brightening profile for a plane-parallel atmosphere |$f(\mu _{\rm em}) = \ln (1 + \mu _{\rm em}^{-1})$| (Haardt 1993). We do not consider an isotropic angular emissivity since for the c-modes which are mostly incompressible perturbations, such an angular profile would result in negligible variability.

The emission angle ϑem = cos −1μem is defined as the spatial angle in the local frame between the emitted photon and the normal to the surface of the accretion disc. The unit normal |$\hat{n}$| is determined using a scalar surface function F(xα) and the projected spatial derivative in the frame comoving with the disc material
\begin{equation} n_\mu \equiv \nabla _\mu F({\boldsymbol x}^{\alpha }) + u_\mu u^\nu \nabla _\nu F({\boldsymbol x}^\alpha ), \qquad \hat{n}_\mu \equiv \frac{n_\mu }{n_\nu n^{\nu }}, \end{equation}
(32)
where uν is the background four-velocity of the disc material. For the unperturbed case, we can use the surface function |$F({\boldsymbol x}) = r\cos \theta = 0$|⁠, which defines the surface of the equatorial plane, giving |$\hat{n}^\nu _o = -(1/r) \theta ^\nu$|⁠, where θν is the unit vector in the θ direction at the equatorial plane. The emission angle is then given as
\begin{equation} \mu _{\rm em} = \cos \vartheta _{\rm em} \equiv \frac{\hat{n}_\mu k^{\mu }}{u_{\nu } k^{\nu }}, \end{equation}
(33)
where all quantities are evaluated at the emission point of the photon at the disc surface.
To extend this to the c-mode we define a vertical perturbation of the surface ξz = ξz(t, r, ϕ). The normal vector is now defined using the surface function
\begin{equation} F({\boldsymbol x}^\alpha ) = r \cos \theta - \xi ^z(t, r, \phi ) = 0. \end{equation}
(34)
Since the velocity of the perturbation is small compared to the background Keplerian flow, we can continue to use the approximation uν ≃ uoν. For simplicity, we will also approximate the surface displacement ξz with the vertical Lagrangian displacement at the disc mid-plane. We also limit our calculations to small perturbations in the limb brightening case such that this approximation does not produce singular values for the angular emissivity.

4 RESULTS

Figs 3 and 4 show constant redshift (g = 1/[1 + z]) contours, for example, systems with black hole spin a = 0.001 and 0.5, respectively. Modulation of the emission at any particular location will modify the observed intensity in that redshift bin. Also shown with dashed lines are the radii corresponding to the IVR for particular modes (n = 0, 1, 2 for a = 0.001 and n = 2 for a = 0.5). The propagation region for each of the corrugation modes is from the ISCO to the IVR. In Figs 6– 8, we calculate the fractional intensity fvar emitted from the variability region r < rIVR for each spectral bin for 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 region is very close to the black hole where gravitational redshift strongly dominates the Doppler shift.

Constant observed frequency ratio (g = 1/[1 + z]) contours for impact parameters $(\hat{\alpha }, \hat{\beta })$ for a disc around a black hole with spin a/M = 0.001 and inclination μo = 0.1, 0.5, 0.7 from left to right. The dotted lines represent the inner vertical resonance, rIVR, for n = 0, 1, 2 perturbations. The variability of the c-modes is confined between the rISCO and rIVR, thus, the spectral variation is confined to the redshift bins which have contours within the mode propagation region.
Figure 3.

Constant observed frequency ratio (g = 1/[1 + z]) contours for impact parameters |$(\hat{\alpha }, \hat{\beta })$| for a disc around a black hole with spin a/M = 0.001 and inclination μo = 0.1, 0.5, 0.7 from left to right. The dotted lines represent the inner vertical resonance, rIVR, for n = 0, 1, 2 perturbations. The variability of the c-modes is confined between the rISCO and rIVR, thus, the spectral variation is confined to the redshift bins which have contours within the mode propagation region.

Open in new tabDownload slide
Constant observed frequency ratio (g = 1/[1 + z]) contours for impact parameters $(\hat{\alpha }, \hat{\beta })$ for a disc around a black hole with spin a/M = 0.5 and inclination μo = 0.1, 0.5, 0.7 from left to right. The dotted lines represent the inner vertical resonance, rIVR, for n = 2 perturbations. In this region, close to rISCO, the gravitational redshift dominates the orbital Doppler boosting. Since the c-mode propagation region is small for high-spin black holes, the redshift range with spectral variability will be correspondingly small.
Figure 4.

Constant observed frequency ratio (g = 1/[1 + z]) contours for impact parameters |$(\hat{\alpha }, \hat{\beta })$| for a disc around a black hole with spin a/M = 0.5 and inclination μo = 0.1, 0.5, 0.7 from left to right. The dotted lines represent the inner vertical resonance, rIVR, for n = 2 perturbations. In this region, close to rISCO, the gravitational redshift dominates the orbital Doppler boosting. Since the c-mode propagation region is small for high-spin black holes, the redshift range with spectral variability will be correspondingly small.

Open in new tabDownload slide

For small amplitude perturbations spectral intensity plots were derived by binning the redshifts for each pixel, with number of bins dependent on the number of pixels subtended by the disc. Example spectra for various phases are shown in Fig. 5, showing the spectral variability for different phases. This spectral variability is highlighted in Figs 9– 13. In these figures, for a given black hole spin parameter a, observer inclination μo = cos θobs, and angular emissivity law we show the expected unperturbed spectrum of the broadened fluorescence line. We also calculate the normalized difference between the broadened line spectrum at a given oscillation phase and the phase-averaged mean spectrum, the spectral variation and plot it as a function of oscillation phase and observed frequency ratio, g, for the simple tilted disc and the n = 0, 1 and 2 c-modes. We normalize the spectral variation here since the amplitude of the perturbations is an arbitrary parameter in our formalism.

The relative intensity (arbitrary scale) of line emission for system with black hole spin a = 0.001, and inclination μo = 0.7, as a function of observed frequency ratio 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 ≃ 0.7M located at r ≃ 6.8M.
Figure 5.

The relative intensity (arbitrary scale) of line emission for system with black hole spin a = 0.001, and inclination μo = 0.7, as a function of observed frequency ratio 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 ≃ 0.7M located at r ≃ 6.8M.

Open in new tabDownload slide

5 DISCUSSION

Using a semi-analytic relativistic raytracing code for the Kerr metric we have calculated the time-dependent line broadening signature of discoseismic corrugation modes in black hole accretion discs. We have shown that detailed spectral timing of the Fe–Kα line variability could demonstrate if such c-modes are the source of LFQPOs observed in accreting black hole systems.

The spectral variability due to c-modes occurs in redshift ranges determined by the propagation region of the mode. In contrast, the simple tilted-precessing disc model shows variability across the entire spectral range of the broadened line which can be compared to Ingram & Done (2012).

The redshift ranges where the variability of the corrugations modes manifests can be easily understood through the g-contours of Figs 3 and 4. In Fig. 3 for a = 0.001 the propagation regions for the c-modes rISCO < r < rIVR, span the range of constant g-contours. This can be directly compared to the spectral range of the variability in Fig. 9. For higher spins, gravitational redshift of the line can begin to dominate the Doppler boosting due to Keplerian motion. The spectral range of the variability is reduced for a = 0.5, as shown in Fig. 12, and drastically reduced for the a = 0.9 (Fig. 13), where the propagation region is much smaller and strongly dominated by gravitational redshift close to the black hole. The spectral range for variability in the a = 0.5 case can be compared directly to the g-contours spanned by the propagation region for the n = 2 mode in Fig. 4.

While in principle the amplitude of the c-modes can be large, we have limited our detailed calculations to small amplitude perturbations as this simplifies the raytracing required to capture the effects on the Fe–Kα emission. If amplitudes of the c-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 qualitatively similar for larger amplitude oscillations. Large amplitude c-modes can result in order unity variation of μem, or even self-shadowing of the oscillation region, and are expected to lead to variations of the order of ∼fvar of the flux in each spectral bin (see Figs 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 the black hole, where rISCO < r < rIVR. For bright sources, such as GRS 1915 + 105, such large variations should be detectable given sufficient spectral (∼0.05 keV) and temporal resolution (∼0.1 s). For high-spin cases, the total fractional flux of the predicted variability decreases significantly, as the variability region shrinks (see Figs 3 and 4). Thus, the c-mode variability of the iron line for high-spin systems would be very difficult to detect. Additionally, we note that c-modes in high-spin systems (a ≳ 0.5) would have frequencies above most observed LFQPOs.

The fraction of the intensity (fvar) that is emitted from the variability region for each observed frequency ratio (g ≡ 1/[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−3, 10−2, 10−1, 0.5, 0.9). For the most redshifted (lowest g) parts of each line profile, the intensity is dominated by emission from the variability region. In this case, we consider angular emissivity to be governed by the limb darkening law, f(μem = (1 + 2.06 μem), with radial emissivity given by ${\cal R}(r) \propto r^{-3}$. For systems with large amplitude perturbations we can expect fractional intensity variation of the order of fvar for a given redshift bin.
Figure 6.

The fraction of the intensity (fvar) that is emitted from the variability region for each observed frequency ratio (g ≡ 1/[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−3, 10−2, 10−1, 0.5, 0.9). For the most redshifted (lowest g) parts of each line profile, the intensity is dominated by emission from the variability region. In this case, we consider angular emissivity to be governed by the limb darkening law, f(μem = (1 + 2.06 μem), with radial emissivity given by |${\cal R}(r) \propto r^{-3}$|⁠. For systems with large amplitude perturbations we can expect fractional intensity variation of the order of fvar for a given redshift bin.

Open in new tabDownload slide
The fraction of the intensity (fvar) that is emitted from the variability region for each observed frequency ratio (g ≡ 1/[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−3, 10−2, 10−1, 0.5, 0.9). For the most redshifted (lowest g) parts of each line profile, the intensity is dominated by emission from the variability region. In this case, we consider angular emissivity to be governed by the limb darkening law, f(μem = (1 + 2.06 μem), with radial emissivity given by ${\cal R}(r) \propto r^{-3}$. For systems with large amplitude perturbations, we can expect fractional intensity variation of the order of fvar for a given g bin.
Figure 7.

The fraction of the intensity (fvar) that is emitted from the variability region for each observed frequency ratio (g ≡ 1/[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−3, 10−2, 10−1, 0.5, 0.9). For the most redshifted (lowest g) parts of each line profile, the intensity is dominated by emission from the variability region. In this case, we consider angular emissivity to be governed by the limb darkening law, f(μem = (1 + 2.06 μem), with radial emissivity given by |${\cal R}(r) \propto r^{-3}$|⁠. For systems with large amplitude perturbations, we can expect fractional intensity variation of the order of fvar for a given g bin.

Open in new tabDownload slide
The fraction of the intensity (fvar) that is emitted from the variability region for each observed frequency ratio (g ≡ 1/[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−3, 10−2, 10−1, 0.5, 0.9). For the most redshifted (lowest g) parts of each line profile, the intensity is dominated by emission from the variability region. In this case, we consider angular emissivity to be governed by the limb darkening law, f(μem = (1 + 2.06 μem), with radial emissivity given by ${\cal R}(r) \propto r^{-3}$. For systems with large amplitude perturbations we can expect fractional intensity variation of the order of fvar for a given g bin.
Figure 8.

The fraction of the intensity (fvar) that is emitted from the variability region for each observed frequency ratio (g ≡ 1/[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−3, 10−2, 10−1, 0.5, 0.9). For the most redshifted (lowest g) parts of each line profile, the intensity is dominated by emission from the variability region. In this case, we consider angular emissivity to be governed by the limb darkening law, f(μem = (1 + 2.06 μem), with radial emissivity given by |${\cal R}(r) \propto r^{-3}$|⁠. For systems with large amplitude perturbations we can expect fractional intensity variation of the order of fvar for a given g bin.

Open in new tabDownload slide

The signature of the radial order (n) of the oscillation is evident primarily in the red and blue edges of the spectral variation. The number of nodes in the oscillation mode is reflected in the number of nodes in the red or blue wings of the spectral variation at a particular phase, as seen in Figs 9– 13. This is particularly prominent for the low-spin cases where the nodes are well spaced (Figs 9– 11). However, such fine signatures in the spectral variation would be difficult to detect unless sufficiently high spectral and timing resolution observations are stacked by oscillation phase.

Upper panels: results assuming limb darkening angular emissivity for a black hole with spin a = 0.001 and inclinations μo = cos θobs = 0.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity.
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.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity.

Open in new tabDownload slide
Upper panels: results assuming limb darkening angular emissivity for a black hole with spin a = 0.01 and inclinations μo = cos θobs = 0.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity.
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.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity.

Open in new tabDownload slide
Upper panels: results assuming limb darkening angular emissivity for a black hole with spin a = 0.1 and inclinations μo = cos θobs = 0.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity.
Figure 11.

Upper panels: results assuming limb darkening angular emissivity for a black hole with spin a = 0.1 and inclinations μo = cos θobs = 0.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity.

Open in new tabDownload slide
Upper panels: results assuming limb darkening angular emissivity for a black hole with spin a = 0.5 and inclinations μo = cos θobs = 0.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity.
Figure 12.

Upper panels: results assuming limb darkening angular emissivity for a black hole with spin a = 0.5 and inclinations μo = cos θobs = 0.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity.

Open in new tabDownload slide
Upper panels: results assuming limb darkening angular emissivity for a black hole with spin a = 0.9 and inclinations μo = cos θobs = 0.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity. The patchiness of the tilted disc spectral variation in this case is due to mainly the resolution and redshift binning and is not physical.
Figure 13.

Upper panels: results assuming limb darkening angular emissivity for a black hole with spin a = 0.9 and inclinations μo = cos θobs = 0.1, 0.5 and 0.7. The first row shows the unperturbed broadened line spectra with the outer edge of the disc rout = 20M. The remaining rows show the spectral variation, the difference in intensity from the average spectra, as a function of oscillation phase and observed frequency ratio for the simple tilted disc, and the n = 0, 1 and 2 c-modes. Lower panels: same as above, but for the limb brightening angular emissivity. The patchiness of the tilted disc spectral variation in this case is due to mainly the resolution and redshift binning and is not physical.

Open in new tabDownload slide

The presence or absence of such signatures in LFQPO-phase stacked Fe–Kα observations can confirm or rule out c-modes as a source of LFQPOs. However, even if variability of the inner disc structure is not the source of the broad-band LFQPOs, observations with high temporal and spectral resolution X-ray instruments, such as on the proposed LOFT or ATHENA missions, would allow iron line probes of inner disc structure variability. The spectral ranges and frequencies over which this variability is seen for c-modes can act as a probe of black hole spin parameter complementary to existing spin constraints from thermal and broadened iron-line observations.

Our quantitative results have relied on raytracing for small amplitude oscillations. Detailed raytracing of large amplitude vertical displacements (ξz ≳ H) was beyond the scope of this paper, however, the scaled spectral variability for large amplitudes should 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 modelling of the oscillation of more realistic discs and their vertical structure should be undertaken to constrain the parameters involved.

DT was supported at Caltech by the Sherman Fairchild Foundation and at McGill through funding from the Canadian Institute for Advanced Research and the Lorne Trottier Chair in Astrophysics and Cosmology. IB conducted this research as part of a Summer Undergraduate Research Fellowship, and was supported by NASA ATP grant no. NNX11AC37G, NSF grant no. PHY-1151197, the David and Lucile Packard Foundation, the Alfred P. Sloan Foundation, Mrs Albert Burford and the 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.

1

This code is available for download at http://www.tapir.caltech.edu/dtsang/qpotrace.html.

2

We calculate the Jacobi elliptic functions and elliptic integrals utilizing standard recurrence relations modified to work with values on part of the complex plane, combined with the Landen's transforms to change the complex arguments.

REFERENCES

Abramowitz
M.
, 
Stegun
I. A.
. , 
Handbook of Mathematical Functions
, 
1965
New York
Dover Publications

Google Scholar

OpenURL Placeholder Text

 
Bardeen
J. M.
, 
Petterson
J. A.
. , 
ApJ
, 
1975
, vol. 
195
 pg. 
L65
 
Crossref
Search ADS
 
Beckwith
K.
, 
Done
C.
. , 
MNRAS
, 
2004
, vol. 
352
 pg. 
353
 
Crossref
Search ADS
 
Čadež
A.
, 
Brajnik
M.
, 
Gomboc
A.
, 
Calvani
M.
, 
Fanton
C.
. , 
A&A
, vol. 
403
 pg. 
29
 
Crossref
Search ADS
 
Cunningham
C. T.
, 
Bardeen
J. M.
. , 
ApJ
, 
1973
, vol. 
183
 pg. 
237
 
Crossref
Search ADS
 
Dexter
J.
, 
Agol
E.
. , 
ApJ
, 
2009
, vol. 
696
 pg. 
1616
 
Crossref
Search ADS
 
Drasco
S.
, 
Hughes
S. A.
. , 
Phys. Rev. D
, 
2004
, vol. 
69
 pg. 
044015
 
Crossref
Search ADS
 
Feroci
M.
, et al. , 
Exp. Astron.
, 
2012
, vol. 
34
 pg. 
415
 
Crossref
Search ADS
 
Fragile
P. C.
, 
Miller
W. A.
, 
Vandernoot
E.
. , 
ApJ
, 
2005
, vol. 
635
 pg. 
157
 
Crossref
Search ADS
 
Fragile
P. C.
, 
Blaes
O. M.
, 
Anninos
P.
, 
Salmonson
J. D.
. , 
ApJ
, 
2007
, vol. 
668
 pg. 
417
 
Crossref
Search ADS
 
Haardt
F.
. , 
ApJ
, 
1993
, vol. 
413
 pg. 
680
 
Crossref
Search ADS
 
Ingram
A.
, 
Done
C.
. , 
MNRAS
, 
2012
, vol. 
427
 pg. 
934
 
Crossref
Search ADS
 
Ipser
J. R.
, 
Lindblom
L.
. , 
ApJ
, 
1992
, vol. 
389
 pg. 
392
 
Crossref
Search ADS
 
Karas
V.
, 
Martocchia
A.
, 
Subr
L.
. , 
PASJ
, 
2001
, vol. 
53
 pg. 
189
 
 
Kato
S.
. , 
PASJ
, 
2001
, vol. 
53
 pg. 
1
 
 
Kato
S.
, 
Fukue
J.
. , 
PASJ
, 
1980
, vol. 
32
 pg. 
377
 
 
Lai
D.
, 
Tsang
D.
. , 
MNRAS
, 
2009
, vol. 
393
 pg. 
979
 
Crossref
Search ADS
 
Laor
A.
. , 
ApJ
, 
1991
, vol. 
376
 pg. 
90
 
Crossref
Search ADS
 
Machida
M.
, 
Matsumoto
R.
. , 
PASJ
, 
2008
, vol. 
60
 pg. 
613
 
 
Miller
J. M.
, 
Homan
J.
. , 
ApJ
, 
2005
, vol. 
618
 pg. 
L107
 
Crossref
Search ADS
 
Misner
C. W.
, 
Thorne
K. S.
, 
Wheeler
J. A.
. , 
Gravitation.
, 
1973
San Francisco
Freeman and Co.

Google Scholar

OpenURL Placeholder Text

 
Noble
S. C.
, 
Krolik
J. H.
, 
Hawley
J. F.
. , 
ApJ
, 
2010
, vol. 
711
 pg. 
959
 
Crossref
Search ADS
 
Novikov
I. D.
, 
Thorne
K. S.
. 
Dewitt
C.
, 
Dewitt
B. S.
. , 
Black Holes (Les Astres Occlus)
, 
1973
New York
Gordon and Breach
pg. 
343
 

Google Scholar

OpenURL Placeholder Text

 
Nowak
M. A.
, 
Wagoner
R. V.
. , 
ApJ
, 
1991
, vol. 
378
 pg. 
656
 
Crossref
Search ADS
 
Nowak
M. A.
, 
Wagoner
R. V.
. , 
ApJ
, 
1992
, vol. 
393
 pg. 
697
 
Crossref
Search ADS
 
Okazaki
A. T.
, 
Kato
S.
, 
Fukue
J.
. , 
PASJ
, 
1987
, vol. 
39
 pg. 
457
 
 
Penna
R. F.
, 
McKinney
J. C.
, 
Narayan
R.
, 
Tchekhovskoy
A.
, 
Shafee
R.
, 
McClintock
J. E.
. , 
MNRAS
, 
2010
, vol. 
408
 pg. 
752
 
Crossref
Search ADS
 
Penna
R. F.
, 
Sä Dowski
A.
, 
McKinney
J. C.
. , 
MNRAS
, 
2012
, vol. 
420
 pg. 
684
 
Crossref
Search ADS
 
Perez
C. A.
, 
Silbergleit
A. S.
, 
Wagoner
R. V.
, 
Lehr
D. E.
. , 
ApJ
, 
1997
, vol. 
476
 pg. 
589
 
Crossref
Search ADS
 
Press
W. H.
, 
Teukolsky
S. A.
, 
Vetterling
W. T.
, 
Flannery
B. P.
. , 
Numerical Recipes in C
, 
1992
Cambridge
Cambridge Univ. Press

Google Scholar

OpenURL Placeholder Text

 
Schnittman
J. D.
, 
Homan
J.
, 
Miller
J. M.
. , 
ApJ
, 
2006
, vol. 
642
 pg. 
420
 
Crossref
Search ADS
 
Silbergleit
A. S.
, 
Wagoner
R. V.
, 
Ortega-Rodríguez
M.
. , 
ApJ
, 
2001
, vol. 
548
 pg. 
335
  
(SWO)
Crossref
Search ADS
 
Sobczak
G. J.
, 
McClintock
J. E.
, 
Remillard
R. A.
, 
Cui
W.
, 
Levine
A. M.
, 
Morgan
E. H.
, 
Orosz
J. A.
, 
Bailyn
C. D.
. , 
ApJ
, 
2000
, vol. 
531
 pg. 
537
 
Crossref
Search ADS
 
Svoboda
J.
, 
Dovčiak
M.
, 
Goosmann
R.
, 
Karas
V.
. , 
A&A
, 
2009
, vol. 
507
 pg. 
1
 
Crossref
Search ADS
 
Svoboda
J.
, 
Dovčiak
M.
, 
Goosmann
R. W.
, 
Jethwa
P.
, 
Karas
V.
, 
Miniutti
G.
, 
Guainazzi
M.
. , 
2012
 
preprint (arXiv:e-prints)
Swank
J. H.
. , 
Nucl. Phys. B Proc. Suppl.
, 
1999
, vol. 
69
 pg. 
12
 
Crossref
Search ADS
 
Tombesi
F.
, 
de Marco
B.
, 
Iwasawa
K.
, 
Cappi
M.
, 
Dadina
M.
, 
Ponti
G.
, 
Miniutti
G.
, 
Palumbo
G. G. C.
. , 
A&A
, 
2007
, vol. 
467
 pg. 
1057
 
Crossref
Search ADS
 
Tsang
D. C.-W.
. , 
PhD thesis
, 
2009
Cornell University

Google Scholar

OpenURL Placeholder Text

 
Tsang
D.
, 
Lai
D.
. , 
MNRAS
, 
2008
, vol. 
387
 pg. 
446
 
Crossref
Search ADS
 
Tsang
D.
, 
Lai
D.
. , 
MNRAS
, 
2009
, vol. 
393
 pg. 
992
 
Crossref
Search ADS
 
Varnière
P.
, 
Tagger
M.
, 
Rodriguez
J.
. , 
A&A
, 
2012
, vol. 
545
 pg. 
A40
 
Crossref
Search ADS
 
Wagoner
R. V.
. , 
Phys. Rep.
, 
1999
, vol. 
311
 pg. 
259
 
Crossref
Search ADS
 

APPENDIX A: RayTracing

In order to calculate various spectrum and light-curve properties we must first construct a simulated image of the black hole and 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 impact parameters α(⊥ to the spin axis projection) and β (∥ to the spin axis projection) disc. Solving for the geodesics, each of these rays can be backtraced 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 G = c = MBH = 1 (e.g. Misner, Thorne & Wheeler 1973)
\begin{eqnarray} \biggl (\frac{{\rm d}t}{{\rm d}\lambda }\biggr ) = \rho ^{-2}\biggl [\frac{r^2 + a^2}{\Delta }[E(r^2 + a^2) - L_z] - a(aE\sin ^2\theta - L_z)\biggr ],\nonumber\\ \end{eqnarray}
(A1)
\begin{eqnarray} \bigg (\frac{{\rm d}r}{{\rm d}\lambda }\biggr )= \rho ^{-2}[(E(r^2 + a^2)-L_za)^2 - \Delta ((L_z - aE)^2 + Q)]^{1/2}, \nonumber\\ \end{eqnarray}
(A2)
\begin{eqnarray} \biggl (\frac{{\rm d}\theta }{{\rm d}\lambda }\biggr ) = \rho ^{-2}[Q - \cos ^2\theta (L_z^2\csc ^2\theta - E^2a^2)]^{1/2}, \end{eqnarray}
(A3)
\begin{eqnarray} \biggl (\frac{{\rm d}\phi }{{\rm d}\lambda }\biggr ) = \rho ^{-2}\biggl [ -aE + L_z \csc ^2\theta + \frac{a}{\Delta }(E(r^2+a^2) - L_za)\biggr ], \nonumber\\ \end{eqnarray}
(A4)
where a is the black hole spin, λ is the affine parameter and Δ = r2 − 2r + a2. E, the photon energy, Lz, the angular momentum 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 authors as 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 & Bardeen (1973) outline a quicker method of calculating the photon trajectories using a Hamilton–Jacobi method. Here we follow a related procedure. Rearranging and switching variables from the affine parameter to a ‘Mino parameter’ (see e.g. Drasco & Hughes 2004) |$\lambda ^{\prime }: \frac{{\rm d}\lambda }{{\rm d}\lambda ^{\prime }} = \rho ^{-2}E$| we get the following coupled first order ODEs for the coordinates as a function of Mino parameter,
\begin{eqnarray} \biggl (\frac{{\rm d}t}{{\rm d}\lambda ^{\prime }}\biggr ) = T(r,\theta ) & \equiv & \biggl [\frac{(r^2 + a^2)^2}{\Delta } - a^2\sin ^2\theta \biggr ]\nonumber\\ &&+\, a l \biggl [1 - \frac{r^2 + a^2}{\Delta }\biggr ], \end{eqnarray}
(A5)
\begin{eqnarray} \bigg (\frac{{\rm d}r}{{\rm d}\lambda ^{\prime }}\biggr )^2 = R(r) & \equiv & (r^2 + a^2 -la)^2 - (r^2 -2r + a^2)\nonumber\\ && \times\,\left[(l-a)^2 +q^2\right], \end{eqnarray}
(A6)
\begin{eqnarray} \biggl (\frac{{\rm d}\theta }{{\rm d}\lambda ^{\prime }}\biggr )^2 = \Theta (\theta ) \equiv q^2 - l^2\cot ^2\theta + a^2 \cos ^2\theta , \end{eqnarray}
(A7)
\begin{eqnarray} \biggl (\frac{{\rm d}\phi }{{\rm d}\lambda ^{\prime }}\biggr ) = \Phi (r,\theta ) \equiv l \csc ^2\theta + a \biggl (\frac{r^2 + a^2}{\Delta } - 1 \biggr ) - \frac{a^2 l}{\Delta }, \end{eqnarray}
(A8)
where l = Lz/E, q2 = Q/E2. The constants of motion l and q2 are related to the impact parameters (α, β) by |$l= -\alpha \sqrt{1-\mu _o^2}$| and |$q^2 = \beta ^2 + \mu _o^2(\alpha ^2 - a^2)$|⁠. By utilizing the ‘Mino parameter’ we avoid code complications involving integrating over multiple turning 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 a semi-analytic fashion using the Jacobi elliptic functions and elliptic integrals. We first solve for the two independent variables, r(λ′) and θ(λ′).

A1 The R equation

We define the roots of the quartic equation
\begin{eqnarray} R(r) = (r^2 + a^2 - al)^2 - (r^2 - 2r + a^2)[(l-a)^2 + q^2]=0,\nonumber\\ \end{eqnarray}
(A9)
 r1, r2, r3 and r4 (see e.g. Čadež et al. 2003), allowing us to rewrite the differential equation for dr/dλ′ as
\begin{eqnarray} \frac{{\rm d}r}{{\rm d}\lambda ^{\prime }} = \sqrt{R(r)} \end{eqnarray}
(A10)
\begin{eqnarray} \hphantom{\frac{{\rm d}r}{{\rm d}\lambda ^{\prime }} }= -\sqrt{(r-r_1)(r-r_2)(r-r_4)(r-r_3)} \end{eqnarray}
(A11)
\begin{eqnarray} \Delta \lambda ^{\prime } = -\int _{r_o}^r \frac{{\rm d}r}{\sqrt{(r-r_1)(r-r_2)(r-r_4)(r-r_3)}} \end{eqnarray}
(A12)
\begin{eqnarray} \hphantom{\Delta \lambda ^{\prime }} &=& \frac{2}{\sqrt{(r_2-r_4)(r_1-r_3)}}{\rm F}\bigg (\sin ^{-1}\sqrt{\frac{(r-r_1)(r_2-r_4)}{(r-r_2)(r_1-r_4)}},\nonumber\\ && \qquad\qquad\qquad\qquad \sqrt{\frac{(r_1-r_4)(r_2-r_3)}{(r_2-r_4)(r_1-r_3)}} \bigg )\bigg |_r^\infty \end{eqnarray}
(A13)
\begin{eqnarray} \hphantom{\Delta \lambda ^{\prime }}&=& \frac{2}{\sqrt{(r_2-r_4)(r_1-r_3)}} {\rm sn}^{-1}\bigg ( \sqrt{\frac{(r-r_1)(r_2-r_4)}{(r-r_2)(r_1-r_4)}},\nonumber\\ && \qquad\qquad\qquad\qquad\sqrt{\frac{(r_1-r_4)(r_2-r_3)}{(r_2-r_4)(r_1-r_3)}} \bigg ) \bigg |_r^\infty , \end{eqnarray}
(A14)
where sn is the Jacobi elliptic function and F is the Jacobi elliptic integral of the first kind.2
With the negative value of the momentum corresponding to backtraced photon decreasing in r. Solving for r(Δλ′) we get
\begin{equation} r(\Delta \lambda ^{\prime }) = \frac{r_1(r_2-r_4)-r_2(r_1-r_4){\rm sn}^2(u, \kappa _r)}{(r_2-r_4) - (r_1-r_4){\rm sn}^2(u, \kappa _r)}, \end{equation}
(A15)
where
\begin{eqnarray} \kappa _r = \sqrt{\frac{(r_1 - r_4)(r_2-r_3)}{(r_2-r_4)(r_1-r_3)}}, \end{eqnarray}
(A16)
\begin{eqnarray} u = \frac{\sqrt{r_2-r_4)(r_1-r_3)}\Delta \lambda ^{\prime }}{2} - u_\infty , \end{eqnarray}
(A17)
\begin{eqnarray} u_\infty = {\rm sn}^{-1}\bigg (\sqrt{\frac{r_2 - r_4}{r_1 - r_4}}, \kappa _r \bigg ) \end{eqnarray}
(A18)
even for complex values of rn.
To calculate the value of Δλ′(r) we must carefully consider any turning points that may be encountered in r. If any root rn is a positive real value greater than the horizon radius then the photon may have a turning point in r. If no turning point is encountered before the emission point then the value of Δλ′(r) is given by
\begin{eqnarray} \Delta \lambda ^{\prime }(r) &=& \frac{2}{\sqrt{(r_2-r_4)(r_1-r_3)}}\biggl [{\rm F}\bigg (\sin ^{-1}\psi _\infty , \kappa _r\bigg )\nonumber\\ && \qquad\quad -\,{\rm F}\bigg (\sin ^{-1}\psi (r), \kappa _r\bigg )\biggr ], \end{eqnarray}
(A19)
where |$\psi (r) = \sqrt{\frac{(r-r_1)(r_2-r_4)}{(r-r_2)(r_1-r_4)}}$| and |$\psi _\infty = \psi (r)|_{r\rightarrow \infty } =\sqrt{\frac{(r_2-r_4)}{(r_1-r_4)}}$|⁠.
If a photon turning point is encountered by the backtrace before reaching the emission point (i.e. Δλ′ > Δλ′(rturn)) then the corresponding value of Δλ′ is given by
\begin{eqnarray} \Delta \lambda ^{\prime }(r) &=& \frac{2}{\sqrt{(r_2-r_4)(r_1-r_3)}}\bigg [{\rm F}\bigg (\sin ^{-1}\psi (r), \kappa _r \bigg )\bigg |_{r_{{\rm turn}}}^\infty \nonumber\\ && \qquad\quad + \, {\rm F}\bigg (\sin ^{-1}\psi (r), \kappa _r \bigg )\bigg |_{r_{{\rm turn}}}^{r_{\rm em}}\bigg ]. \end{eqnarray}
(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 2θ such that
\begin{eqnarray} \frac{{\rm d}\theta }{{\rm d}\lambda ^{\prime }} = \pm \sqrt{\frac{-a^2z^2 - z[q^2 +l^2 -a^2] + q^2}{1- z}} \end{eqnarray}
(A21)
\begin{eqnarray} \hphantom{\frac{{\rm d}\theta }{{\rm d}\lambda ^{\prime }}}= \pm \sqrt{\frac{a^2(z_+ - z)(z - z_-)}{1-z}}, \end{eqnarray}
(A22)
where |$z_\pm = -\frac{q^2 + l^2 -a^2}{2a^2} \pm \sqrt{\frac{(q^2 + l^2 -a^2)^2}{4a^4} + \frac{q^2}{a^2}}$|⁠. If we take χ: z = z+cos 2χ we have
\begin{equation} \frac{{\rm d}\chi }{{\rm d}\theta } = \pm \sqrt{\frac{1-z}{(z_+ - z)}}, \end{equation}
(A23)
which gives
\begin{eqnarray} \frac{{\rm d}\chi }{{\rm d}\lambda ^{\prime }} = \frac{{\rm d}\chi }{{\rm d}\theta }\frac{{\rm d}\theta }{{\rm d}\lambda ^{\prime }} = \sqrt{a^2(z - z_-)} \end{eqnarray}
(A24)
\begin{eqnarray} \hphantom{\frac{{\rm d}\theta }{{\rm d}\lambda ^{\prime }}}= \sqrt{a^2(z_+ \cos ^2\chi -z_-)}. \end{eqnarray}
(A25)
This gives the integral
\begin{eqnarray} \lambda ^{\prime }(\chi ) - \lambda ^{\prime }(\chi _o) = \int _{\chi _o}^\chi \frac{{\rm d}\chi }{\sqrt{a^2(z_+\cos ^2\chi - z_-)}} \end{eqnarray}
(A26)
\begin{eqnarray} \quad = \frac{1}{a\sqrt{z+ - z_-}}\biggl [{\rm F}\biggl (\chi , \sqrt{\frac{z_+}{z_+-z_-}}\biggr )- {\rm F}\biggl (\chi _o, \sqrt{\frac{z_+}{z_+-z_-}}\biggr ) \biggr ],\nonumber\\ \end{eqnarray}
(A27)
where F(ψ, k) is the elliptic integral of the first kind. (In Abramowitz and Stegun notation this is F = F(ψ|m), where m = k2.)
Inverting this equation and solving for θ we finally obtain
\begin{equation} \theta (\lambda ^{\prime }) = \cos ^{-1}(\sqrt{z_+} {\rm cn}(a\sqrt{z_+ - z_-}\lambda ^{\prime } + u_{\theta _o}, \kappa _\theta )), \end{equation}
(A28)
where
\begin{eqnarray} u_{\theta _o} = {\rm sgn}(\beta ){\rm cn}^{-1}(\cos (\theta _o)/\sqrt{z_+}, \kappa _\theta ), \end{eqnarray}
(A29)
\begin{eqnarray} \kappa _\theta = \sqrt{\frac{z_+}{z_+ - z_-}}, \end{eqnarray}
(A30)
and cn(u, κ) is the Jacobi elliptic function which has inverse cn− 1(u, κ) = F(cos − 1(u), κ).
In order to calculate the Mino parameter corresponding to a particular value of θ we see
\begin{eqnarray} \Delta \lambda ^{\prime }(\theta ) &=& \frac{1}{a\sqrt{z_+ - z_-}}\biggl [{\rm F}\biggl (\chi , \sqrt{\frac{z_+}{z_+-z_-}}\biggr )\nonumber\\ && \qquad\quad -\,{\rm F}\biggl (\chi _o, \sqrt{\frac{z_+}{z_+-z_-}}\biggr )\biggr ], \end{eqnarray}
(A31)
where
\begin{eqnarray} \chi = \cos ^{-1}\biggl (\frac{\cos \theta }{\sqrt{z_+}}\biggr ), \end{eqnarray}
(A32)
\begin{eqnarray} \chi _o = {\rm sgn}(\beta )\cos ^{-1}\biggl (\frac{\cos \theta _o}{\sqrt{z_+}}\biggr ), \end{eqnarray}
(A33)
thus, we can solve for λ′ corresponding to intersection of the ray with simple fixed θem disc configurations. The flat disc model of θem = π/2 is of particular interest.

A3 The Φ equation

The ϕ(λ′) differential equation is more complicated than the equations for the first two spatial coordinates; however, the differential equation can be solved by breaking the integration into two parts, integration over θ and integration over r.

First rewriting the Φ(r, θ) equation we see
\begin{equation} \frac{{\rm d}\phi }{{\rm d}\lambda ^{\prime }} = \frac{l}{1-\cos ^2\theta } + a\frac{2r + la}{r^2 - 2r + a^2}. \end{equation}
(A34)
Integrating the first term we see
\begin{eqnarray} \Delta \phi _1 = \int _{\lambda ^{\prime }_o}^{\lambda _{\rm em}} \frac{l {\rm d}\lambda ^{\prime }}{1 - \cos ^2 \theta } \end{eqnarray}
(A35)
\begin{eqnarray} \hphantom{\Delta \phi _1}= \int _{\lambda ^{\prime }_o}^{\lambda ^{\prime }_{\rm em}}\frac{l {\rm d}\lambda ^{\prime }}{1 - z_+{\rm cn}^2(a\sqrt{z_+-z_-}\lambda ^{\prime } + u_{\theta _o}, k_\theta )} \end{eqnarray}
(A36)
\begin{eqnarray} \hphantom{\Delta \phi _1}= \int _{\lambda ^{\prime }_o}^{\lambda _{\rm em}} \frac{l{\rm d}\lambda ^{\prime }}{1 - z_+[1 - {\rm sn}^2(u, k_\theta )]} \end{eqnarray}
(A37)
\begin{eqnarray} \hphantom{\Delta \phi _1}= \frac{l}{(1-z_+)}\int _{\lambda ^{\prime }_o}^{\lambda ^{\prime }_{\rm em}}\frac{{\rm d}\lambda ^{\prime }}{1 + \frac{z_+}{1-z_+}{\rm sn}^2(u, k_\theta )} \end{eqnarray}
(A38)
\begin{eqnarray} \hphantom{\Delta \phi _1}= \frac{l/a}{(1-z_+)\sqrt{z_+-z_-}}\int _{u_o}^{u_{\rm em}}\frac{{\rm d}u}{1 + \frac{z_+}{1-z_+}{\rm sn}^2(u, k_\theta )} \end{eqnarray}
(A39)
\begin{eqnarray} \hphantom{\Delta \phi _1}= \frac{l/a}{(1-z_+)\sqrt{z_+-z_-}}\Pi \bigg (\chi , \frac{z_+}{1-z_+}, \kappa _\theta \bigg )\bigg |_{\chi _o}^{\chi _{\rm em}}, \end{eqnarray}
(A40)
where |$u = a\sqrt{z_+-z_-}\lambda ^{\prime } + u_o$|⁠, χ = am(u, κθ) is the Jacobi amplitude of u and Π(ψ, n, κ) is the elliptic integral of the third kind. 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:
\begin{eqnarray} \Delta \phi _2 = a\int _{\lambda ^{\prime }_o}^{\lambda ^{\prime }_{\rm em}} \frac{(2r + la){\rm d}\lambda ^{\prime }}{r^2 - 2r + a^2} \end{eqnarray}
(A41)
\begin{eqnarray} \hphantom{\Delta \phi _2 }= a\int ^{r_{\rm em}}_{\infty }\frac{2r + la}{(r-r_+)(r-r_-)}\frac{{\rm d}\lambda ^{\prime }}{{\rm d}r}{\rm d}r \end{eqnarray}
(A42)
\begin{eqnarray} \hphantom{\Delta \phi _2 }&=& a\bigg (\frac{2r_+ + al}{r_+-r_-}\bigg )\int _{r_{\rm em}}^{\infty }\frac{1}{(r-r_+)}\nonumber\\ && \times\,\frac{{\rm d}r}{\sqrt{(r-r_1)(r-r_2)(r-r_3)(r-r_4)}}-\, a\bigg (\frac{2r_- + al}{r_+-r_-}\bigg )\nonumber \\ && \times\,\int _{r_{\rm em}}^{\infty }\frac{1}{(r-r_-)}\frac{{\rm d}r}{\sqrt{(r-r_1)(r-r_2)(r-r_3)(r-r_4)}} \nonumber\\ \end{eqnarray}
(A43)
\begin{eqnarray} \hphantom{\Delta \phi _2 }&=& a\bigg [\frac{2r_2 + al}{r_2^2-2r_2 + a^2}\bigg ]\Delta \lambda ^{\prime }+ \frac{2a}{\sqrt{(r_1-r_3)(r_2-r_4)}}\bigg (\frac{r_1-r_2}{r_+-r_-}\bigg )\nonumber \\ && \times \bigg [\frac{2r_-+al}{(r_1-r_-)(r_2-r_-)}\Pi _- - \frac{2r_++al}{(r_1-r_+)(r_2-r_+)}\Pi _+\bigg ]_{r_{\rm em}}^{\infty },\nonumber\\ \end{eqnarray}
(A44)
where
\begin{eqnarray} r_\pm = 1 \pm \sqrt{1-a^2} \end{eqnarray}
(A45)
\begin{eqnarray} \Pi _\pm = \Pi (\psi (r), n_\pm , \kappa _r) \end{eqnarray}
(A46)
\begin{eqnarray} n_\pm = \psi ^2(r)\bigg |_{r=r_\pm }. \end{eqnarray}
(A47)
Thus, we have
\begin{equation} \Delta \phi (\Delta \lambda ^{\prime }) = \Delta \phi _1 + \Delta \phi _2. \end{equation}
(A48)
Note that when a turning point is encountered in r, the Δϕ2 must be evaluated with the appropriate sign change as for Δλ′(r).

In general, the parameters of the elliptical integrals of the third kind will be complex. This limits evaluation of Δϕ2 using recurrence relations to a particular domain of complex parameter space. Evaluations outside of this domain can be performed using 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
\begin{equation} \frac{{\rm d}t}{{\rm d}\lambda ^{\prime }} = - a^2\sin ^2\theta + \frac{(r^2 + a^2)^2 + 2alr}{r^2 - 2r + a^2}. \end{equation}
(A49)
The first term can be integrated in a relatively straightforward fashion as
\begin{eqnarray} \Delta t_1 = \int a^2(1 + \cos ^2\theta ) \mathrm{d}\lambda ^{\prime } \end{eqnarray}
(A50)
\begin{eqnarray} \hphantom{\Delta t_1}= a^2\Delta \lambda ^{\prime } - \int a^2 z_+ {\rm cn}^2(u_\theta , \kappa _\theta ))\frac{{\rm d}u_\theta }{a\sqrt{z_+-z_-}} \end{eqnarray}
(A51)
\begin{eqnarray} \hphantom{\Delta t_1}= a^2 \Delta \lambda ^{\prime } - \frac{az_+}{\sqrt{z_+-z_-}}\frac{1}{\kappa _\theta ^2}\bigg [{\rm E}({\rm am}(u_\theta ), \kappa _\theta ) - (1-\kappa _\theta ^2)u_\theta \bigg ]_{u_o}^{u_{\rm em}} \nonumber\\ \end{eqnarray}
(A52)
\begin{eqnarray} \hphantom{\Delta t_1}= a^2\bigg [a + \frac{z_+(1-\kappa _\theta ^2)}{\kappa _\theta ^2}\bigg ]\Delta \lambda ^{\prime } - \frac{az_+}{\kappa _\theta ^2\sqrt{z_+-z_-}}{\rm E}(\chi , \kappa _\theta )\bigg |_{\chi _o}^{\chi _{\rm em}}. \nonumber\\ \end{eqnarray}
(A53)

For the second term, |$\Delta t_2\equiv \frac{(r^2 + a^2)^2 + 2alr}{r^2 - 2r + a^2}$|⁠, the analytic solution can be determined as a very long combination of elliptical integrals 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 are necessarily divergent. Though one could simply evaluate Δt2 only up to a large arbitrary value of r, it is better to instead subtract off the same infinite constant for each ray, as our interest is limited to the elapsed observer time difference between different rays, by considering the Kerr time.

Remembering that the time in Kerr coordinates is given by a transformation
\begin{equation} {\rm d}t_{\rm Kerr} = {\rm d}t_{\rm BL} + \bigg (\frac{r^2 + a^2}{r^2-2r+a^2}\bigg ) \mathrm{d}r, \end{equation}
(A54)
we can then subtract off the same constant for each ray by taking
\begin{eqnarray} \Delta t^{\prime } = \Delta t_{\rm Kerr} - \int _{r_{\rm em}}^{r_{\rm bitrary}}\frac{r^2 + a^2}{r^2-2r + a^2} \mathrm{d}r \end{eqnarray}
(A55)
\begin{eqnarray} \hphantom{\Delta t^{\prime }}= \Delta t_{\rm BL} + \int _\infty ^{r_{\rm em}} \frac{r^2 + a^2}{r^2 - 2r + a^2}dr - \int _{r_{\rm em}}^{r_{\rm bitrary}}\frac{r^2 + a^2}{r^2-2r + a^2} \mathrm{d}r \nonumber\\ \end{eqnarray}
(A56)
\begin{eqnarray} \hphantom{\Delta t^{\prime }}= \Delta t_1 + \Delta t_2 - \int _{r_{\rm bitrary}}^\infty \frac{r^2 + a^2}{r^2 -2r + a^2}\mathrm{d}r \end{eqnarray}
(A57)
\begin{eqnarray} \hphantom{\Delta t^{\prime }}&=& \Delta t_1 + \int _{r_{\rm em}}^{r_{\rm bitrary}} \frac{(r^2 + a^2)^2 + 2alr}{r^2 - 2r + a^2}\frac{{\rm d}r}{\sqrt{R(r)}} \nonumber \\ &&+\, \int _{r_{\rm bitrary}}^\infty \bigg [\frac{(r^2 + a^2)^2 + 2alr}{r^2 - 2r + a^2}\frac{1}{\sqrt{R(r)}} - \frac{r^2 + a^2}{r^2 - 2r + a^2}\bigg ]\mathrm{d}r,\nonumber\\ \end{eqnarray}
(A58)
where rbitrary is some r beyond the disc range.

These integrations are best performed with numerical quadrature as evaluation of the elliptic integrals and Jacobi elliptic functions prove less efficient than the numerical integration. In addition, the |$1/\sqrt{R(r)}$| term makes it difficult to remove the divergent components of the analytic solution so it may be evaluated numerically.

© 2013 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
Issue Section:
Article
Download all slides