ALGORITHMS AND PROGRAMS FOR STRONG GRAVITATIONAL LENSING IN KERR SPACE-TIME INCLUDING POLARIZATION

Accretion disks are one of the most spectacular phenomena in modern astrophysics. The large luminosity of active galactic nuclei (AGNs) and quasars (as large as 10⁴⁷ erg s⁻¹) is believed to be the result of gas accreted by central supermassive black holes. Up to ∼10% of the accreting mass can be converted to radiated energy during the accretion process (much higher than in nuclear fusion processes, about ∼0.5%). Besides emission in the infrared, optical, and UV bands, most AGNs emit X-rays. Unlike optical emission, which can be explained by standard accretion disk theory (Shakura & Sunyaev 1973), the physical mechanism of AGN X-ray emission remains an enigma, despite decades of research (Krolik 2008). According to standard theory, the temperature of an AGN accretion disk is not hot enough to emit X-rays. The X-ray emission is thought to be generated through inverse Compton scattering of disk photons with hot electrons in the so-called X-ray "corona." However, both the physics and geometry of this mysterious X-ray corona are not well understood. Observations have revealed that photons of different frequencies correspond to emission regions of different sizes: the higher the frequency, the smaller the emission size, and the closer it is to the central black hole. Very recently, observations in quasar microlensing (Kochanek 2004; Blackburne et al. 2006, 2014, 2015; Morgan et al. 2008, 2012; Chartas et al. 2009; Dai et al. 2010; Chen et al. 2011, 2012; Mosquera et al. 2013; MacLeod et al. 2015) have constrained the AGN X-ray emission size to be of the order of 10 gravitational radii (${{r}_{g}}\equiv GM/{{c}^{2}}$, where M is the mass of the central black hole). Within distances of a few ${{r}_{g}},$ Einstein's general relativistic gravity theory plays a fundamental role in understanding accreting phenomena. In particular, the X-ray emission is strongly lensed by the central black hole before it escapes the gravity field and arrives at a distant observer. It is well known that photons follow null geodesics in curved spacetimes. The spacetime around a rotating black hole is described by the Kerr metric (Kerr 1963). The strong field environment of AGN X-ray emission invalidates the linear approximation as made in standard gravitational lensing theory (Schneider et al. 1992). Consequently, numerical integration of the geodesic equations is necessary for many problems involving the Kerr metric. Another well-known fact is that a photon's polarization vector is parallel transported along the photon's geodesic and its plane of polarization is rotated as the photon passes the rotating black hole, i.e., the so-called gravitational Faraday rotation (Balazs 1958; Plebanski 1960; Ishihara et al. 1988; Agol 1997; Agol & Krolik 2000; Frolov & Shoom 2012; Yoo 2012). Since X-rays are emitted in regions very close to the black hole, the observed X-ray polarization will be influenced by the strong gravity field. Following the pioneering work of Stark, Connors, and Piran (Connors & Stark 1977; Stark & Connors 1977; Connors et al. 1980), X-ray polarization has now become a very promising probe of important parameters such as black hole spin and the inclination angle of the accretion disk (Agol 1997; Agol & Krolik 2000; Dovčiak et al. 2004; Li et al. 2009; Schnittman et al. 2013). It is therefore essential to include the effects of strong gravity when studying AGN X-ray emission.

Numerical ray-tracing in the Kerr metric is not new (Cunningham 1975; Luminet 1979; Rauch & Blandford 1994; Bromley et al. 1997; Beckwith & Done 2004; Fuerst & Wu 2004; Schnittman & Bertschinger 2004; Broderick 2006; Cadeau et al. 2007; Dolence et al. 2009; Psaltis & Johannsen 2012; Chan et al. 2013; Chen et al. 2013a, 2013b; Schnittman & Krolik 2013). However, public codes are relatively rare. Exceptions are Dexter & Agol (2009; a FORTRAN code, fast but mathematically convolved; see also Yang & Wang 2013), Vincent et al. (2011; a C++ code, not parallelized), Kuchelmeister et al. (2012), and Chan et al. (2013; both are GPU-based ray-tracing code solving second-order geodesic equations using CUDA). Polarization and gravitational Faraday rotation have been studied extensively in the literature (Agol 1997; Agol & Krolik 2000; Schnittman & Krolik 2009); however, there seems to be no public code with a focus on X-ray polarizations. Furthermore, there seems to be no ray-tracing code written in interpreted high-level languages such as MATLAB or Python. We present simple algorithms for computing the effect of strong lensing of AGN accretion disks by a Kerr black hole, including effects on polarization. Our algorithms are realized in the high-level languages MATLAB and Python, and they are fully parallelized and can be easily realized in other languages such as FORTRAN and C++.

The plan of this paper is as follows. An introduction to the Kerr metric, the tetrad formalism, the strong lensing formalism, and the polarization formalism is given in Section 2. The graphic user interface (GUI) and numerical algorithms are discussed in detail in Section 3. In Section 4, we give an example of the use of our MATLAB/Python program KERTAP (KErr Ray-Tracing And Polarization). In Section 5, we conclude.

Since the gas accreted by black holes carries angular momentum, the black holes powering the accretion processes are expected to possess angular momentum themselves. The proper metric describing the spacetime around a rotating black hole was discovered by Kerr in 1963 (Kerr 1963) and was named after him. We introduce the Kerr metric in this section and discuss the other formalism needed in the strong lensing and polarization analysis of AGN X-ray emission.

2.1. Kerr Metric

In Boyer–Lindquist coordinates, the Kerr metric can be written as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} d{{s}^{2}} & = & -{{\alpha }^{2}}d{{t}^{2}}+\frac{{{{\Sigma }}^{2}}}{{{\rho }^{2}}}{{{\rm sin} }^{2}}\theta {{(d\phi -\omega dt)}^{2}} \\ {} & {} & +\frac{{{\rho }^{2}}}{{\Delta }}d{{r}^{2}}+{{\rho }^{2}}d{{\theta }^{2}}, \\ \end{array}\end{eqnarray} \tag{ 1 }$$

where $(t,r,\theta ,\phi )$ are the four independent coordinate variables, and we have defined

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {{\rho }^{2}} & \equiv & {{r}^{2}}+{{a}^{2}}{{{\rm cos} }^{2}}\theta \\ {\Delta } & \equiv & {{r}^{2}}-2Mr+{{a}^{2}} \\ {{{\Sigma }}^{2}} & \equiv & {{\rho }^{2}}\left[ ({{r}^{2}}+{{a}^{2}})+\frac{2Mr}{{{\rho }^{2}}}{{a}^{2}}\;{{{\rm sin} }^{2}}\theta \right] \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\alpha }^{2}} & \equiv & \frac{{{\rho }^{2}}{\Delta }}{{{{\Sigma }}^{2}}} \\ \omega & \equiv & \frac{2aMr}{{{{\Sigma }}^{2}}} \\ \end{array}\end{eqnarray} \tag{ 2 }$$

with the constants M and a being the mass and angular momentum (i.e., J/M) of the Kerr black hole (we have taken $c=G=1$). The geodesic motion of a photon (its 4 momentum satisfies ${{p}^{a}}{{p}_{a}}=0$, index a runs from 1 to 4, with 1 being the time component) can be described by a group of eight first-order Hamilton equations (Arnold 1978):

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{d{{x}^{a}}}{d\xi } & = & \frac{\partial {\mathcal{H}}}{\partial {{p}_{a}}} \\ \frac{d{{p}_{a}}}{d\xi } & = & -\frac{\partial {\mathcal{H}}}{\partial {{x}^{a}}}, \\ \end{array}\end{eqnarray} \tag{ 3 }$$

where ξ is the affine parameter and ${\mathcal{H}}({{x}^{a}},{{p}_{b}})\equiv \frac{1}{2}{{g}^{ab}}(x){{p}_{a}}{{p}_{b}}$ is the Hamiltonian. The number of equations to integrate is reduced to six since p_t and ${{p}_{\phi }}$ are motion constants (the Hamiltonian ${\mathcal{H}}$ does not depend on t or ϕ). Another motion constant was discovered by Carter (1968):

$$\begin{eqnarray}&&{\mathcal{C}}=p_{\theta }^{2}+{{{\rm cos} }^{2}}\theta \left( \frac{p_{\phi }^{2}}{{{{\rm sin} }^{2}}\theta }-{{a}^{2}}p_{t}^{2} \right).\end{eqnarray} \tag{ 4 }$$

The number of ordinary differential equations (ODEs) used to solve the geodesic equations could have been reduced by one using this constant; however, we use Carter's constant ${\mathcal{C}}$ as an independent check of our ray-tracing code. A Kerr black hole has two horizons (outer and inner) defined as ${{r}_{\pm }}=M\pm \sqrt{{{M}^{2}}-{{a}^{2}}}$. The maximum angular momentum of a Kerr black hole was found in Thorne (1974) to be ${{a}_{{\rm max} }}=0.998.$ We stop the ray-tracing when the geodesic is within $0.1\%$ of the outer horizon ${{r}_{+}}.$ For accretion disks (normally assumed in the equatorial plane, $\theta =\pi /2$), an important quantity is the so-called innermost-stable-circular-orbit (ISCO) in the equatorial plane (Chandrasekhar 1983). For example, ${{r}_{{\rm ISCO}}}=6\;{{r}_{g}}$ and $1.24\;{{r}_{g}}$, respectively, for a Schwarzschild black hole (a = 0) and an extreme Kerr black hole. We use ${{r}_{{\rm ISCO}}}$ as the inner cutoff of the accretion disk in this paper.

2.2. Tetrad Formalism

In curved spacetime, the frequency (or energy) of a photon depends on the spacetime point where the frequency is measured, and also on the 4 velocity u^a of the observer measuring the frequency (${{u}^{a}}{{u}_{a}}=-1$). If p_a is the 4 momentum of the photon, then the observed photon frequency $\nu =-{{u}^{a}}{{p}_{a}}.$ Let ${{\nu }_{{\rm o}}}$ and ${{\nu }_{{\rm e}}}$ be the photon's frequency measured by a distant observer and an observer co-moving with the light source (emitter). We define the redshift factor

$$\begin{eqnarray}&&g\equiv \frac{{{\nu }_{{\rm o}}}}{{{\nu }_{e}}}=\frac{{{(-{{u}^{a}}{{p}_{a}})}_{{\rm o}}}}{{{(-{{u}^{a}}{{p}_{a}})}_{{\rm e}}}},\end{eqnarray} \tag{ 5 }$$

where $u_{{\rm o}}^{a}$ and $u_{{\rm e}}^{a}$ are the respective 4 velocities of the distant observer and the source.

For the distant observer, we choose ${{u}^{a}}=\frac{1}{\alpha }(1,0,0,\omega ),$ the so-called zero-angular-momentum observer (ZAMO; Bardeen et al. 1972), where ω and α are defined in Equation (2). A photon's momentum can be specified using the orthonormal tetrad of ZAMO instead of the coordinate basis, i.e.,

$$\begin{eqnarray}&&\left( \begin{array}{ccccccccccccccc} {{{\hat{p}}}^{t}} \\ {{{\hat{p}}}^{r}} \\ {{{\hat{p}}}^{\theta }} \\ {{{\hat{p}}}^{\phi }} \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} \alpha & 0 & 0 & 0 \\ 0 & \frac{\rho }{\sqrt{{\Delta }}} & 0 & 0 \\ 0 & 0 & \rho & 0 \\ -\frac{{\Sigma }\;{\rm sin} \;\theta }{\rho }\omega & 0 & 0 & \frac{{\Sigma }\;{\rm sin} \;\theta }{\rho } \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} {{p}^{t}} \\ {{p}^{r}} \\ {{p}^{\theta }} \\ {{p}^{\phi }} \\ \end{array} \right).\end{eqnarray} \tag{ 6 }$$

We will use Equation (6) to compute the polarization measured by a distant observer (see Section 2.5). The 4 velocity of the source is model-dependent. For accretion disks, we assume a Keplerian flow for the accreting gas (Bardeen et al. 1972), i.e.,

$$\begin{eqnarray}&&u_{K}^{a}=\gamma (1,0,0,{{{\Omega }}_{K}}),\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{{\Omega }}_{K}} & \equiv & \frac{{{M}^{1/2}}}{{{r}^{3/2}}+a{{M}^{1/2}}} \\ \gamma & \equiv & {{(-{{g}_{11}}-2{{g}_{14}}{{{\Omega }}_{K}}-{{g}_{44}}{\Omega }_{K}^{2})}^{-1/2}}. \\ \end{array}\end{eqnarray} \tag{ 8 }$$

We similarly define an orthonormal tetrad in the comoving frame of the Keplerian flow, and to change between tetrad and coordinate components we use

$$\begin{eqnarray}&&\left( \begin{array}{ccccccccccccccc} {{p}^{t}} \\ {{p}^{r}} \\ {{p}^{\theta }} \\ {{p}^{\phi }} \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} \gamma & 0 & 0 & \tau \lambda \\ 0 & \frac{\sqrt{{\Delta }}}{\rho } & 0 & 0 \\ 0 & 0 & -\frac{1}{\rho } & 0 \\ \gamma {{{\Omega }}_{K}} & 0 & 0 & \tau \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} {{{\hat{p}}}^{1}} \\ {{{\hat{p}}}^{2}} \\ {{{\hat{p}}}^{3}} \\ {{{\hat{p}}}^{4}} \\ \end{array} \right),\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \lambda & \equiv & -\frac{{{g}_{14}}+{{g}_{44}}{{{\Omega }}_{K}}}{{{g}_{11}}+{{g}_{14}}{{{\Omega }}_{K}}} \\ \tau & \equiv & {{({{g}_{11}}{{\lambda }^{2}}+2{{g}_{14}}\lambda +{{g}_{44}})}^{-1/2}}. \\ \end{array}\end{eqnarray} \tag{ 10 }$$

The above equations are useful since both the intensity and the polarization profile of the source emission are often given in the local rest frame of the source. In Section 2.5, we use Equations (9) to initialize the polarization propagator at the end of backward ray-tracing.

2.3. Strong Lensing by Kerr Black Holes

X-ray emission from an accretion disk is strongly lensed by the central Kerr black hole as it leaves the source to arrive at a distant observer. A lensed accretion disk appears very different from an unlensed one. The frequency of the source emission is gravitationally and Doppler blue- or redshifted depending on the position of the emitter, the 4 velocity of the source, and the distance and inclination angle of the observer. Consequently, the intensity profile of a lensed accretion disk can be very different from an unlensed one, and the observed flux can also be different from the unlensed case.

By definition, the monochromatic flux ${{F}_{{{\nu }_{{\rm o}}}}},$ measured by an observer O very far from the black hole, is (Mihalas 1978)

$$\begin{eqnarray}&&{{F}_{{{\nu }_{{\rm o}}}}}\equiv \oint {{I}_{{{\nu }_{{\rm o}}}}}\mu d{{{\Omega }}_{{\rm o}}}=\int {{I}_{{{\nu }_{{\rm o}}}}}d{{{\Omega }}_{{\rm o}}},\end{eqnarray} \tag{ 11 }$$

where $d{{{\Omega }}_{{\rm o}}}$ is the differential solid angle measured at the observer, and μ is the angle cosine between the light ray and the normal of the detector window. We can safely drop the factor μ for sources at cosmological redshifts because these sources are observed with tiny solid angles. If we ignore Kerr strong lensing and Doppler shifts, photons travel along straight lines and there is no frequency redshift. Consequently, ${{I}_{\nu }}$ is a conserved quantity along the light path,

$$\begin{eqnarray}&&F_{{{\nu }_{{\rm o}}}}^{{\rm unlensed}}\equiv \int I_{{{\nu }_{{\rm o}}}}^{{\rm obs}}({\boldsymbol{\hat{n}}} )d{{{\Omega }}_{{\rm o}}}=\frac{{\rm cos} \;\theta }{r_{{\rm obs}}^{2}}\int I_{{{\nu }_{{\rm o}}}}^{{\rm source}}({\boldsymbol{\hat{n}}} )d{\mathcal{A}},\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{\hat{n}}}$ is the 3D photon direction at the observer, $d{\mathcal{A}}$ is the differential area element of the accretion disk, and we have changed the solid angle integration at the observer into a 2D surface integral over the accretion disk. On the other hand, for Kerr lensing, we have

$$\begin{eqnarray}&&F_{{{\nu }_{{\rm o}}}}^{{\rm lensed}}\equiv \int {{g}^{3}}({\boldsymbol{\hat{n}}} )I_{{{\nu }_{{\rm e}}}}^{{\rm source}}({{x}^{a}}({\boldsymbol{\hat{n}}} ),{{p}_{a}}({\boldsymbol{\hat{n}}} ))d{{{\Omega }}_{{\rm o}}},\end{eqnarray} \tag{ 13 }$$

where ${{\nu }_{{\rm e}}}$ is the source frequency, g is the redshift factor of Equation (5), and $({{x}^{a}},{{p}_{a}})$ is the photon's position and momentum at the emitter found by backward ray-tracing (see Section 2.4). To obtain Equation (13), we simply used the fact that ${{I}_{\nu }}/{{\nu }^{3}}$ is conserved along each geodesic. The strong lensing magnification of the specific flux is then

$$\begin{eqnarray}&&\mu \equiv \frac{F_{{{\nu }_{{\rm o}}}}^{{\rm lensed}}}{F_{{{\nu }_{{\rm o}}}}^{{\rm unlensed}}}.\end{eqnarray} \tag{ 14 }$$

2.4. Backward Ray-tracing

2.5. Polarization and Gravitational Faraday Rotation

Because the polarization at the source is known and needs to be computed at the observer, most current polarization propagators are based on forward ray-tracing (Schnittman & Krolik 2009). Our polarization integrator is based on backward ray-tracing. We assume the polarization and limb darkening of the source emission are given as in the classical work of Chandrasekhar (1960) where the radiation transfer equation with polarization was solved assuming a scattering-dominated semi-infinite atmosphere. Other initial polarization profiles such as an inverse Compton scattering-dominated optically thin atmosphere can also be implemented.

At the end of the backward ray-tracing of each ray, we know the landing point on the accretion disk of a geodesic and its 4 momentum p^a. First, we compute the photon's 4 momentum components ${{\hat{p}}^{a}}$ in the comoving frame. We then compute the angle cosine, $\mu ,$ between the photon's momentum and the disk normal measured by the comoving observer. The angular-dependent part of the intensity profile, $w(\mu )$ (see Equation (21) in Section 4), and the degree of polarization, ${{\delta }_{{\rm source}}}(\mu )$, can be computed by interpolating Table XXIV of Chandrasekhar (1960). We then compute the photon's polarization vector ${\boldsymbol{E}}$ in the local rest frame using the fact that ${\boldsymbol{E}}$ is parallel to the disk plane and normal to the photon's momentum vector ($E\cdot p=0$), i.e.,

$$\begin{eqnarray}&&{{\hat{E}}^{a}}=\frac{1}{\sqrt{{{({{{\hat{p}}}^{2}})}^{2}}+{{({{{\hat{p}}}^{4}})}^{2}}}}(0,-{{\hat{p}}^{4}},0,{{\hat{p}}^{2}}).\end{eqnarray} \tag{ 15 }$$

Next, we use Equation (9) to compute the coordinate components ${{E}^{a}}$ of the polarization vector ${\boldsymbol{E}} .$ We then solve the parallel propagation equation of the polarization vector forward along the photon geodesic toward the distant observer using

$$\begin{eqnarray}&&\frac{d{{E}^{a}}}{d\xi }=-{\Gamma }_{\,bc}^{a}{{p}^{b}}{{E}^{c}},\end{eqnarray} \tag{ 16 }$$

where ${\Gamma }_{bc}^{a}$ are the Christoffel symbols. Since the geodesic equations have already been solved backward from the observer to the source, and because we have dropped those rays that missed the accretion disk, the integral from the source forward to the observer is straightforward and fast.

At the end of the forward polarization propagating process, we are back to the observer with rotated polarization vector ${{{\boldsymbol{E}} }_{{\rm obs}}}$. Here, we first compute the tetrad components of ${{{\boldsymbol{E}} }_{{\rm obs}}}$ using Equation (6). We then work out the Stokes parameters corresponding to this ray,

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} 2\chi & = & {\rm arctan} \left( -\frac{{{{\hat{E}}}^{\theta }}}{{{{\hat{E}}}_{\phi }}} \right) \\ Q & = & \delta {{I}_{{{\nu }_{o}}}}{\rm cos} 2\chi \\ U & = & \delta {{I}_{{{\nu }_{o}}}}{\rm sin} 2\chi , \\ \end{array}\end{eqnarray} \tag{ 17 }$$

where χ is the polarization angle, $\delta ={{\delta }_{{\rm source}}}$ is the degree of polarization conserved along a photon geodesic, and ${{I}_{{{\nu }_{o}}}}$ is the observed specific intensity at frequency ${{\nu }_{o}},$ which is related to the source intensity by ${{I}_{{{\nu }_{o}}}}={{g}^{3}}{{I}_{{{\nu }_{e}}}}$. To compute the integrated degree of polarization and the polarization angle, we add $Q,$ $U,$ and ${{I}_{{{\nu }_{o}}}}$ from each image pixel, obtaining ${{Q}_{{\rm total}}}$, ${{U}_{{\rm total}}},$ and $I_{{{\nu }_{o}}}^{{\rm total}}.$ Finally, we use Equation (17) again to compute ${{\chi }_{{\rm total}}}$ and ${{\delta }_{{\rm total}}}.$ Detailed algorithms are given in the next section.

3.1. Graphic User Interface

The MATLAB version of KERTAP contains a GUI (see Figure 1). It asks the user for the black hole parameters, the resolution desired, and the number of parallel workers for parallel computing. The GUI contains a status panel indicating the status of the code: waiting for input, running, or finished. If the code is run successfully, the GUI will generate the redshifted image of the accretion disk, and output important parameters such as the flux magnification and the degree and angle of the observed polarization. The performance of the code, such as the wall clock time used for ray-tracing, will also be recorded.

Zoom In Zoom Out Reset image size

3.2. Initialization Algorithm For Backward Ray-tracing

Given an observer distance ${{r}_{{\rm obs}}},$ and the dimensions of the image rectangle (orthogonal to the line of sight at the source redshift), ${{x}_{{\rm size}}}$ and ${{y}_{{\rm size}}},$ the angular size of the cone focusing at the observer can be inferred. Given the resolutions in the x and y directions, the pixel size can be computed as $d{{x}_{{\rm grid}}}={{x}_{{\rm size}}}/{{n}_{{\rm x}}}$ and $d{{y}_{{\rm grid}}}={{y}_{{\rm size}}}/{{n}_{{\rm y}}}.$ For a pixel labeled by ($ix,iy$), the photon's starting 8D phase-space coordinates are computed in Algorithm 1.

Algorithm 1.  Initialization of one pixel ($ix,iy$) of the grid. ${{\hat{p}}^{r,\theta ,\phi }}$ are elements of the 3D photon momentum vector with respect to the ZAMO. The function p3_to_p4() will normalize it assuming ${{\hat{p}}^{t}}=1,$ and then convert the ZAMO tetrad components to covariant coordinate components. We choose ${{t}_{{\rm obs}}}=0$ and ${{\phi }_{{\rm obs}}}=0$ without a loss of generality.

${{X}_{0}}=(0,{{r}_{{\rm obs}}},{{\theta }_{{\rm obs}}},0)$	$\triangleright$ observer's 4 coordinates
${{\hat{p}}^{r}}={{r}_{{\rm obs}}}$	$\triangleright$ must be positive
${{\hat{p}}^{\phi }}=d{{x}_{{\rm grid}}}\times ix$
${{\hat{p}}^{\theta }}=d{{y}_{{\rm grid}}}\times iy$
$p4=p3\_to\_p4({{\hat{p}}^{r}},{{\hat{p}}^{\theta }},{{\hat{p}}^{\phi }})$	$\triangleright$ p4 is covariant
${{Y}_{0}}=({{X}_{0}},p4)$	$\triangleright$ 8D starting point for RT

Download table as:  ASCIITypeset image

3.3. Ray-tracing Using MATLAB ODE45

For each ray arriving at the observer, we solve the geodesic equations using MATLAB ODE solver ode45. The syntax of ode45 is

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} [\xi ,Y,{{\xi }_{e}},{{Y}_{e}},{{I}_{e}}] \\ \quad =ode45(Y^{\prime} ,{{\xi }_{{\rm range}}},{{Y}_{0}},{\rm options},\ldots ), \\ \end{array}\end{eqnarray} \tag{ 18 }$$

where $Y^{\prime} (\xi ,Y)$ is an 8D column vector containing the right-hand side of Equation (3), Y₀ is the 8D initial condition for ode45, and ${{\xi }_{{\rm range}}}$ is the range of the independent variable (affine parameter) where we want to solve the ODEs. A subtle point of backward ray-tracing is that we do not shoot rays toward the accretion disk (say, take ${{\hat{p}}^{r}}\lt 0$), instead we must have ${{\hat{p}}^{r}}\gt 0$ (toward the observer). This point is important because the Kerr metric is not symmetric with respect to the reflection of the time coordinate (time reversal will change the direction of the black hole angular momentum). The backward ray-tracing is realized by properly choosing ${{\xi }_{{\rm range}}},$ i.e., we define ${{\xi }_{{\rm range}}}=[{{\xi }_{0}},{{\xi }_{1}}]$, where $Y({{\xi }_{0}})={{Y}_{0}},$ and $Y({{\xi }_{1}})$ corresponds to the light emitting event (${{\xi }_{1}}\lt {{\xi }_{0}}$, and ${{t}_{1}}\lt {{t}_{0}}$). In Equation (18), options is a structure containing optional parameters that change the default integration properties, such as the relative and absolute error tolerance, and events to stop the integration. For example, suppose we want to stop the ray-tracer when the geodesic is very close to the event horizon, say, $r\lt (1+\delta ){{r}_{+}}$ for some small δ, we define a black hole hitting event as

$$\begin{eqnarray}&&{\rm fla}{{{\rm g}}_{{\rm BH}}}=(r\leqslant (1+\delta ){{r}_{+}})-1.\end{eqnarray} \tag{ 19 }$$

Since we want to know when the geodesic hits the equatorial plane, we define another event

$$\begin{eqnarray}&&{\rm fla}{{{\rm g}}_{{\rm Equator}}}={\rm cos} \;\theta .\end{eqnarray} \tag{ 20 }$$

Then, ode45 will record the point where ${\rm fla}{{{\rm g}}_{{\rm BH}}}=0$ or ${\rm fla}{{{\rm g}}_{{\rm Equator}}}=0$ and stop the integration if the user requests to (as we did). These two events are defined in disk_events.m. In the left-hand side of Equation (18), the index I_e tells you which event happens, $({{\xi }_{e}},{{Y}_{e}})$ stores the point where the event happens, and $(\xi ,Y)$ stores the results of the integration up to the event point (i.e., the photon geodesic). Based on ode45, the ray-tracing part of KERTAP is easily realized in Algorithm 2.

Algorithm 2.  Ray-tracing Algorithm. MATLAB ode45 is Based on Runge–Kutta45.

for $ix=-nx:nx$	do
for $iy=-ny:ny$	do
Initialize pixel ($ix,iy$) giving Y0	$\triangleright$ See Algorithm 1
$[\xi ,Y,\ldots ]=ode45(Y^{\prime} (\xi ),{{Y}_{0}},{\rm options},\ldots )$
if ray goes to infinity	then
goes to next ray
els if hit on the horizon
${{n}_{{\rm hole}}}={{n}_{{\rm hole}}}+1$
go to next ray
els if hit on the equatorial plane	then
if hit on the accretion disk	then
${{n}_{{\rm disk}}}={{n}_{{\rm disk}}}+1$
compute redshift, write $({{x}^{a}},{{p}_{b}})$, ...
go to next ray
end if
end if
end for
end for

Download table as:  ASCIITypeset image

3.4. Polarization Propagator

At the end of the backward ray-tracing of each ray, the polarization propagator is initialized and propagated forward to the observer, see Algorithm 3.

Algorithm 3.  Polarization Propagating Algorithm. In the following, polarization() computes the polarization degree ${{\delta }_{{\rm source}}}$ and polarization vector E^a_source on the accretion disk using ${{({{x}^{a}},{{p}^{a}})}_{{\rm source}}}$ obtained from backward ray-tracing. ${{\delta }_{{\rm source}}}$ is interpolated from Table XXIV of Chandrasekhar (1960). The polarization vector E^a is propagated forward from the source to the observer using Equation (16). polar_obs() computes the polarization angle χ at the observer. ${{\delta }_{{\rm obs}}}={{\delta }_{{\rm source}}}$ since the degree of polarization is conserved along geodesics.

for $ix=-nx:nx$	do
for $iy=-ny:ny$	do
$[\xi ,Y,{{Y}_{e}}]=ode45(Y^{\prime} ,\ldots )$	$\triangleright$ Algorithm 2
if a good shot	then
$[{{\delta }_{{\rm source}}},{{E}_{{\rm source}}}]={\rm polarization}({{Y}_{e}},M,a)$
${{E}_{{\rm curr}}}={{E}_{{\rm source}}}$
${{n}_{{\rm point}}}={\rm length}(\xi$)
for $j=0,{{n}_{{\rm point}}}-2$ do	$\triangleright$ ${{E}_{{\rm source}}}\to {{E}_{{\rm obs}}}$
$d\xi ={{\xi }_{j+1}}-{{\xi }_{j}}$
for $a=1:4$	do
$d{{E}^{a}}=-\sum _{b,c=1}^{4}\left( {\Gamma }_{\,bc}^{a}{{p}^{b}}E_{{\rm curr}}^{c} \right)d\xi$
$E_{{\rm next}}^{a}=E_{{\rm curr}}^{a}+d{{E}^{a}}$
end for
${{E}_{{\rm curr}}}={{E}_{{\rm next}}}$
end for
${{E}_{{\rm obs}}}={{E}_{{\rm curr}}}$
$\chi ={\rm polar}\_{\rm obs}({{E}_{{\rm obs}}})$
end if
end for
end for

Download table as:  ASCIITypeset image

3.5. Post Processing

The results of ray-tracing are stored in different data structures. For example, Image_G and Image_Y store the redshifts (g = 0 for a missed shot) and 8D phase-space coordinates $({{x}^{a}},{{p}_{b}})$ of the target points on the accretion disk for all successful shots. Image_P stores the polarization data $(w,\delta ,\chi )$, where w is the limb darkening factor from Table XXIV of Chandrasekhar (1960). The strong lensing amplifications of the observed flux and image area are computed in Algorithm 4. The observed degree and angle of polarization are computed using Algorithm 5. Thanks to the data visualization tools intrinsic to MATLAB and Python, the image of the lensed disk can be easily generated using functions such as imagesc() for MATLAB and imshow() for Python, and the gravitational Faraday rotated polarization vector field can be conveniently visualized using the function quiver(), which is similarly defined for the two languages. The redshifted image of an accretion disk will be plotted in the GUI after a successful run, and another two image windows (one for the same redshift image, the other for the intensity profile and polarization visualization) will pop up for the user to edit or save the figures in their desired formats (eps, jpg, etc.).

Algorithm 4.  Strong lensing algorithm. dA and $d{\Omega }$ are the differential image area and solid angle for one pixel. The unlensed image area and flux are easy to compute, and are assumed to be known.

${{F}^{{\rm lensed}}}=0$; ${{A}^{{\rm lensed}}}=0$	$\triangleright$ Initialization
for $ix=-nx:nx$	do
for $iy=-ny:ny$	do
if a good shot	then
${{A}_{{\rm lensed}}}={{A}_{{\rm lensed}}}+1$
$g={\rm Image}\_{\rm G}(ix,iy)$	$\triangleright$ read the redshift g
$Y={\rm Image}\_{\rm Y}(ix,iy)$	$\triangleright$ read 8D $({{x}^{a}},{{p}_{b}})$
$dF={{g}^{3}}{{I}_{{{\nu }_{o}}/g}}(Y)$	$\triangleright$ ${{\nu }_{e}}={{\nu }_{o}}/g$
${{F}^{{\rm lensed}}}={{F}^{{\rm lensed}}}+dF$
end if
end for
end for
${{A}^{{\rm lensed}}}={{A}^{{\rm lensed}}}\times dA$
${{F}^{{\rm lensed}}}={{F}^{{\rm lensed}}}\times d{\Omega }$
${{\mu }_{{\rm area}}}={{A}^{{\rm lensed}}}/{{A}^{{\rm unlensed}}}$	$\triangleright$ Area amplification
${{\mu }_{{\rm flux}}}={{F}^{{\rm lensed}}}/{{F}^{{\rm unlensed}}}$	$\triangleright$ Flux amplification

Download table as:  ASCIITypeset images: 1 2

Algorithm 5.  Polarization synthesizing algorithm. ${{\delta }_{{\rm total}}}$ and ${{\chi }_{{\rm total}}}$ are the integrated degree and angle of polarization.

${{I}_{{\rm total}}}=0;{{Q}_{{\rm total}}}=0;{{U}_{{\rm total}}}=0$	$\triangleright$ Initialization
for $ix=-nx:nx$	do
for $iy=-ny:ny$	do
if IsImagePoint	then
$g={\rm Image}\_{\rm G}(ix,iy)$	$\triangleright$ read the redshift g
$r={\rm Image}\_{\rm Y}(ix,iy,2)$	$\triangleright$ read r coordinate
$[w,\delta ,\chi ]={\rm Image}\_{\rm P}(ix,iy)$
$dI=w{{g}^{{\Gamma }+2}}/{{r}^{n}}$	$\triangleright$ See Equation (23)
$dQ=\delta \times dI\times {\rm cos} (2\chi )$
$dU=\delta \times dI\times {\rm sin} (2\chi )$
${{Q}_{{\rm total}}}={{Q}_{{\rm total}}}+dQ$
${{U}_{{\rm total}}}={{U}_{{\rm total}}}+dU$
${{I}_{{\rm total}}}={{I}_{{\rm total}}}+dI$
end if
end for
end for
${{\chi }_{{\rm total}}}=0.5{{{\rm tan} }^{-1}}({{U}_{{\rm total}}}/{{Q}_{{\rm total}}})$
${{\delta }_{{\rm total}}}={{U}_{{\rm total}}}/({{I}_{{\rm total}}}{\rm sin} (2{{\chi }_{{\rm total}}}))$

Download table as:  ASCIITypeset image

3.6. Python Version of KERTAP

We assume that the X-ray corona is a thin layer immediately above the AGN accretion disk, and moves with Keplerian flows (see Equations (7) and (8)). We take ${{r}_{{\rm disk}}}=20\;{{r}_{g}}$ according to the constraints obtained from quasar microlensing observations (Kochanek 2004; Blackburne et al. 2006, 2014, 2015; Morgan et al. 2008, 2012; Chartas et al. 2009; Dai et al. 2010; Mosquera et al. 2013; MacLeod et al. 2015), and we use ${{r}_{{\rm ISCO}}}$ as the inner cutoff. The image window for backward ray-tracing is $50\;{{r}_{g}}\times 50\;{{r}_{g}}.$ X-ray coronae of this type are motivated by the so-called sandwich model (Haardt & Maraschi 1991, 1993) and are often used in the literature (Ruszkowski & Fabian 2000). For a corona geometry such as a sphere or a disk above the accretion disk, see Chen et al. (2013a). Because quasar X-ray emission is observed to follow a power law, we assume the following form for the source intensity profile,

$$\begin{eqnarray}&&{{I}_{\nu }}(\nu ,\mu ,r)\propto \frac{1}{{{r}^{n}}}\frac{w(\mu )}{{{\nu }^{{\Gamma }-1}}},\end{eqnarray} \tag{ 21 }$$

where μ is the cosine between the photon's 4 momentum and the upward disk normal measured by the comoving observer, and $w(\mu )$ is the angular-dependence of the intensity profile (we have assumed azimuthal symmetry for simplicity). To be more specific, we take $w(\mu )$ from Chandrasekhar (1960; classical results from solving the radiation transfer equation with polarization assuming a scattering-dominated semi-infinite atmosphere). The parameter Γ is referred to as the photon index, and n specifies the radial steepness of the profile. We take Γ = 2.0 (e.g., Chen et al. 2012), and n = 3 (radially steep) or 0 (radially flat). We test for two black hole spins, a = 0.998 (Thorne 1974) or 0 (Schwarzschild black hole). We take the observer distance as ${{r}_{{\rm obs}}}={{10}^{6}}{{r}_{g}}$ and test for different inclination angles (from nearly face on to nearly edge on).

For the assumed simple geometry and intensity profile, the unlensed specific flux at frequency ${{\nu }_{o}}$ can be analytically integrated using Equation (12)

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} F_{{{\nu }_{o}}}^{{\rm unlensed}} & = & \frac{2\pi {{\mu }_{{\rm obs}}}w({{\mu }_{{\rm obs}}})}{r_{{\rm obs}}^{2}}\frac{1}{\nu _{o}^{{\Gamma }-1}}\int _{{{r}_{{\rm ISCO}}}}^{{{r}_{{\rm disk}}}}\frac{dr}{{{r}^{n-1}}} \\ {} & = & \frac{2\pi {{\mu }_{{\rm obs}}}w({{\mu }_{{\rm obs}}})}{(n-2)r_{{\rm obs}}^{2}}\frac{1}{\nu _{o}^{{\Gamma }-1}}\left[ \frac{1}{r_{{\rm ISCO}}^{n-2}}-\frac{1}{r_{{\rm disk}}^{n-2}} \right], \\ \end{array}\end{eqnarray} \tag{ 22 }$$

where ${{\mu }_{{\rm obs}}}={\rm cos} ({{\theta }_{{\rm obs}}}).$ We have dropped the unimportant constant in the definition of ${{I}_{\nu }}.$ For n = 2, the $(n-2)$ factor in the denominator is absent and the expression within the square brackets is replaced by ${\rm log} \;{{r}_{{\rm disk}}}/{{r}_{{\rm ISCO}}}$. From Equation (13), the lensed flux is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} F_{{{\nu }_{o}}}^{{\rm lensed}} & = & \frac{1}{\nu _{o}^{{\Gamma }-1}}\int \frac{{{g}^{{\Gamma }+2}}w(\mu )}{{{r}^{n}}}d{{{\Omega }}_{o}} \\ {} & = & \frac{1}{r_{{\rm obs}}^{2}}\frac{1}{\nu _{o}^{{\Gamma }-1}}\left( \mathop{\sum }\limits_{i,j}\frac{g_{i,j}^{{\Gamma }+2}w({{\mu }_{i,j}})}{r_{i,j}^{n}} \right)d{{A}_{{\rm grid}}}, \\ \end{array}\end{eqnarray} \tag{ 23 }$$

where the summation is over the accretion disk image, and $d{{A}_{{\rm grid}}}=d{{x}_{{\rm grid}}}\times d{{y}_{{\rm grid}}}$ is the area of one image pixel. Note that $w(\mu )$ varies from pixel to pixel because of the general relativistic light bending and the aberration effect. The magnification $\mu =F_{{{\nu }_{o}}}^{{\rm lensed}}/F_{{{\nu }_{o}}}^{{\rm unlensed}}$ can be computed using Equations (22) and (23) immediately (note that the term $1/(r_{{\rm obs}}^{2}\nu _{o}^{{\Gamma }-1})$ cancels out).

In Figure 2, we plot the lensed intensity profile and distribution of polarization for the X-ray disk observed at inclination angle 75°. Figure 3 is similar to Figure 2, but is for spin a = 0.5 and we have included up to fourth-order Kerr-lensing images. Kerr strong lensing changes both the shape and the intensity profile of the disk significantly. The intensity is strongly concentrated in a small region (in dark red) near the black hole and to the left of the accretion disk (where the rotating source is approaching the observer). As for polarization, the degree of polarization to the left of the black hole is on average smaller than the classical value (∼4.6% for θ = 75°; Chandrasekhar 1960) because the angle between the photon's 4 momentum and the disk's normal measured by the comoving observer is smaller than the inclination angle. Therefore, $\delta (\mu )$ is smaller, see the black lines in the second row of Figure 4 for the inclination dependence of δ. The contrary holds for the right-hand side of the accretion disk. On the other hand, the direction of polarization can change significantly because of the gravitational Faraday rotation. This effect is more significant for source regions closer to the black hole where the observed polarization angle χ can be positive or negative (depending on the actual source location) instead of horizontally aligned as in Chandrasekhar (1960).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The strong lensing magnification of the image area and specific flux, and the integrated (or averaged) polarization degree and angle are plotted in Figure 4. We did not consider higher-order images here since the accretion disk is believed by many to be optically thick. We plot the inclination dependence of ${{\mu }_{{\rm area}}},$ ${{\mu }_{{\rm flux}}},\delta ,$ and χ for spin a = 0.998 and $0,$ and for a steep and flat radial profile (n = 3 or 0). The magnification increases with inclination angles and is more significant for steeper profiles and for larger spins (smaller ${{r}_{{\rm ISCO}}}$), both of which amount to more concentrated emission closer to the black hole. The fact that for small inclination angles (nearly face on), the observed degree of polarization is greater than the classical result is mainly caused by gravitational light bending ($\delta =0$ for the face on case, see Chandrasekhar 1960). For moderate to high inclination angles, the emission is strongly focused in a small region with a low degree of polarization, see Figure 2. Consequently, special and general relativity effects tend to reduce the observed degree of polarization, confirming earlier results (Connors et al. 1980).

Another interesting point is that for a fixed emission size (e.g., ${{r}_{{\rm disk}}}=20\;{{r}_{g}}$ as used in this paper) with the intensity profile Equation (21), the effects on polarization, i.e., the polarization profile, the integrated degree and angle of polarization, are achromatic. This is because both the inputs of the polarization propagator, i.e., $\delta (\mu )$ and $w(\mu )$ from Chandrasekhar (1960), and gravitational Faraday rotation are wavelength independent, and the 1/$\nu _{o}^{{\Gamma }-1}$ factor cancels out when averaging over the image of the disk when computing ${{\chi }_{{\rm total}}}$ and ${{\delta }_{{\rm total}}}$, see Equations (17), (22), and (23). This apparent discrepancy with Schnittman & Krolik (2009) is caused by the fact that Schnittman & Krolik (2009) have assumed thermal emission for X-rays emitted by the accretion disk around a stellar mass black hole, and therefore the X-ray emission of higher frequencies are emitted by regions closer to the black hole (the effective temperature of the disk is higher closer to the black hole), whereas we have assumed the same power-law $\propto {{\nu }^{-({\Gamma }-1)}}$ at each radius for X-ray emission by AGNs. If we assume that hard X-rays are emitted by a region smaller than that of the soft X-rays and closer to the black hole as suggested by Chen et al. (2011), say, ${{r}_{{\rm disk}}}=10\;{{r}_{g}},$ then, even under the same power-law assumption, the integrated degree and angle of polarization should differ from current results, i.e., the results will become chromatic.

The accuracy of KERTAP is very good. We have used Carter's constant as an independent check of accuracy. At the end of the long ray-tracing (${{r}_{{\rm obs}}}={{10}^{6}}{{r}_{g}}$), the relative error of Carter's constant is never worse than $\sim {{10}^{-7}}$ for all cases tried by the authors. We also checked the accuracy of KERTAP using the norm of the photon's 4 momentum vector and the polarization vector E^a; see Figure 5 for the results of a typical ray. The results of the MATLAB code are highly consistent with the Python code, e.g., the norm of the difference between the redshift images generated using the Python and MATLAB code is of the order of $\sim {{10}^{-6}}$, and there were no rays that hit the accretion disk or entered the horizon in one case that did not hit the disk or enter the horizon in the other. A convergence test is shown in Table 1 and Figure 6. The convergence of KERTAP is pretty good, e.g., the amplification to image area and specific flux, and the polarization degree and angle all converge at a grid of moderate size 500 × 500. In Figure 7, we test the scalability of KERTAP. Both the MATLAB (parallelized using MCDS) and Python (parallelized using mpi4py) codes scale roughly linearly with the number of CPUs. Consequently, the performance of KERTAP can be significantly improved by parallel computing. However, the performance improvement of the MATLAB code through parallel computing was restricted by the number of MDCS licenses available in a user's MATLAB cluster, whereas the Python code is only restricted by the size of the cluster (number of CPUs available for the MPI job). It is important to note that since MATLAB2014b, the Parallel Computing Toolbox allows up to 512 parallel workers (this number is not constrained by the number of MDCS licenses available). Consequently, significant speedup can still be achieved by users without a MATLAB cluster but with a single node with many cores.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Convergence and Performance of KERTAP. We Run KERTAP on a Single Computer with Eight MATLAB Parallel Workers

Resolution	${{\mu }_{{\rm area}}}/{{\mu }_{{\rm flux}}}$	${{g}_{{\rm min} }}/{{g}_{{\rm max} }}$	δ	χ (deg)	Wall time (hr)	Speed (ray s−1)
$20\times 20$	1.563/1.742	0.257/1.280	0.0230	−8.473	0.004	27.33
$50\times 50$	1.568/1.587	0.051/1.334	0.0167	−6.728	0.011	66.15
$100\times 100$	1.567/1.712	0.039/1.356	0.0159	−8.017	0.034	83.65
$200\times 200$	1.563/1.723	0.037/1.357	0.0158	−8.052	0.120	93.07
$500\times 500$	1.563/1.726	0.038/1.360	0.0158	−8.005	0.713	97.75
$1000\times 1000$	1.563/1.726	0.036/1.360	0.0158	−7.989	2.881	96.62

Note. Accretion disk size ${{r}_{{\rm disk}}}=20\;{{r}_{g}}.$ The size of the image window is $50\;{{r}_{g}}\times 50\;{{r}_{g}}.$ Black hole spin a=0.998. Photon index Γ = 2.0 and radial profile n = 3. Inclination angle $\theta =75{}^\circ .$ Observer distance ${{r}_{{\rm obs}}}={{10}^{6}}\;{{r}_{g}}$.

Download table as:  ASCIITypeset image

We have developed KERTAP, the first ray-tracing code using the interpreted languages MATLAB and Python, for studying strong gravitational lensing in Kerr spacetime including polarization.³ The key ingredients of KERTAP are a GUI (for the MATLAB version), backward ray-tracing realized by a fifth-order Runge–Kutta method, a polarization propagator and synthesizer, and the formalism and algorithms for computing strong lensing effects of a Kerr black hole centered on an accretion disk. The backward ray-tracing part of KERTAP most resembles that of Schnittman & Bertschinger (2004) and Vincent et al. (2011). Vincent et al. (2011) focuses on radiation transfer while this paper focuses on the observable effects of Kerr strong lensing and polarization. The independence of the geodesic equations from the polarization propagation equations makes it possible to solve the geodesic equations (backward) first, then propagate the polarization vector forward. The idea of backward ray-tracing followed by forward polarization propagating is new, to the best of our knowledge. Polarization and gravitational Faraday rotation have been studied extensively in the literature (e.g., Agol 1997; Agol & Krolik 2000; Schnittman & Krolik 2009). However, public code explicitly dealing with polarization is not yet available. Published papers investigating the X-ray polarization of accretion disks are basically all doing forward ray-tracing (with or without Monte-Carlo radiation transfer; see, e.g., Schnittman & Krolik 2009). Since X-rays are emitted by sources very close to the central black hole and moving with highly relativistic flows, in order to compute the observable property measured by a distant observer at a given inclination angle, at each point on the accretion disk and for each direction, many forward rays have to be traced (because no one knows which ray will go where). This increases the forward ray-tracing time significantly, and the results are usually noisier than those from backward ray-tracing (e.g., compare Figure 1 of Schnittman & Krolik 2009 with Figure 2 of this paper).

The idea of performing very computationally intensive jobs such as ray-tracing in Kerr spacetime using interpreted languages such as MATLAB or Python might sound odd at first. It is hard to expect a dynamic language such as Python to compete in performance with FORTRAN or C++.⁴ In fact, the single core (MATLAB or pure Python) version of KERTAP performs poorly compared to code written in FORTRAN and C++, and to GPU-accelerated CUDA code (see, e.g., Dexter & Agol 2009; Chan et al. 2013). However, this does not mean that interpreted languages such as Python or MATLAB will/should not be used in cpu intensive astronomical computations such as Kerr ray-tracing. The performance of these languages has been improved significantly over the past decade, in particular, through the development of application programming interfaces (APIs) with low-level languages such as C/FORTRAN (e.g., the CPython API), multi-cpu parallel computing (mpi4py for Python; MDCS for MATLAB), and by supporting GPUs. For example, both Python and MATLAB support GPU programming, e.g., PyCUDA⁵ for Python. The performance of the pure Python version of KERTAP was improved by two orders of magnitude using Cython. The disadvantage in speed is compensated by parallel computing. In MATLAB, this is achieved by using the Parallel Computing Toolbox on local machines with multiple cores, or by using the Distributed Computing Server over a cluster. For Python, the Cython version of KERTAP supports full implementation of MPI. For example, in this paper, we have shown that it is feasible to study strong lensing of X-ray emission of accretion disks around a Kerr black hole using MATLAB/Python. The strong lensing analysis of an X-ray emitting disk of radius $20\;{{r}_{g}}$ can be done with high accuracy in ∼40 minutes by running the MATLAB version of KERTAP on a desktop computer with multiple cores (≤1 minute using the Cython code). The computing time can be further reduced on a high-performance computing cluster given the massively parallel nature of ray-tracing and the good scalability of KERTAP (see Figure 7). The MATLAB GUI of KERTAP allows a ray-tracing job to be created on a local machine (e.g., the user's laptop), submitted to a remote high-performance computing cluster with the results sent back to the user upon a successful run, the most desired form of interactive high-performance computing.⁶ If the reader desires, the strong lensing and polarization algorithms can be easily realized in FORTRAN or C++.

The advantage of coding using high-level languages such as Python and MATLAB is multifold. First, the object-oriented nature of MATLAB makes it very easy to build GUIs. For example, the GUI as shown in Figure 1 is written using the MATLAB GUI building tool "guide." Second, it is much easier to write code using MATLAB or Python than FORTRAN or C/C++. Programming in MATLAB/Python greatly shortens the whole cycle of a scientific project, i.e., from prototyping and debugging the code, to production runs, to visualization of the scientific results and generating high quality plots. For example, this whole paper can be done using a single language (either MATLAB or Python). Our code is very short. The kernel of the MATLAB version of KERTAP is only about 600 lines. The kernel of the Python version of KERTAP is about 200 lines longer (including the Runge–Kutta solver). The reader should experience no difficulty in understanding our code after reading this paper. Third, since MATLAB was designed to be a data analysis and visualization tool (and thanks to free Python packages such as "matplotlib"), the post-processing of ray-tracing data using MATLAB or Python is almost trivial. For example, the piece of code generating the redshifted image of the accretion disk (see Figure 1) is only a few lines. Users do not have to write separate codes visualizing their results using software such as IDL or Mathematica. Our code would be very useful for scientists new to the Kerr geometry or in need of polarization calculations, or to those people who need a second code to compare with their own. KERTAP is the first public code explicitly dealing with strong lensing of X-ray polarization. For less complicated computing jobs, our MATLAB/Python code can be used directly or after some modification. For complicated computing jobs, KERTAP might still be useful if the user's institute has a parallel computing cluster since KERTAP is fully parallelized. We conclude that it is possible to perform complex numerical-relativity-related computations using interpreted languages such as MATLAB and Python.

NSF AST-0707704 and support for Program number HST-GO-12948.04-A was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. X.D. acknowledges the financial support from the NASA ADAP programs NNX11AD09G, NNX15AF04G, and NSF grant AST-1413056. We thank the anonymous referee for a careful review of this work. We thank Jeremy D. Schnittman for comments. B.C. thanks Peter Beerli and Paul van der Mark for discussions on parallel computing using Python and Shiyun Tang for help in generating one plot. The performance analysis was conducted on the HPC cluster at the Research Computing Center at the Florida State University.