Accurate closed-form trajectories of light around a Kerr black hole using asymptotic approximants

Light deflection in curved spacetimes is one of the earliest predictions of Einstein’s general theory of relativity, and one of the best understood aspects of the theory. The null geodesics describing photon trajectories have been investigated for a wide variety of physical configurations in a number of limits. Mere years after general relativity was developed, the weak-field properties of light deflection around the sun were used to test the theory by Eddington and others during the eclipse of 1919. The limit where photons approach the innermost circular orbit (ICO¹¹1This applies for massless particles and occurs at different radii than the innermost stable circular orbit [ISCO] that pertains to massive bodies., see figure 1), referred to as the strong-field limit, has also been explored for decades, tracing back to work by Hagihara in the 1930’s [1], and further studies by Darwin [2] in the following decades. After the initial construction of the Kerr metric describing spinning black holes [3], many of the early results on null geodesics in these spacetimes were derived by Carter, beginning in the 1960’s [4]. We refer readers to Chandrasekhar’s work on the subject for a thorough review on geodesics in black hole spacetimes [5].

In our previous paper [6] (hereafter referred to as Paper 1), we examined the deflection of photons traveling in the equatorial plane of a Kerr (spinning) black hole. Building off the work of Iyer and collaborators in particular [7, 8, 9], our work generated approximate expressions for the total light deflection, $\alpha$, as a function of the blackhole (henceforth abbreviated as BH) spin and impact parameter, $b$ (see figure 1). To do so, we utilized asymptotic approximants, i.e., closed-form expressions that connect asymptotic expansions obtained in different limits. These approximants yielded highly accurate expressions for light bending. Asymptotic approximants had previously been used to study thermodynamic phase behavior [10, 11, 12] and the solution of nonlinear boundary value problems [13].

In Paper 1, we considered the total light deflection angle only, which describes the case where a photon begins at infinite (or quasi-infinite) separations from the black hole, passes by, and returns to infinite separations. While gravitational lensing may be examined via the total bending angle, studies involving electromagnetic emission and absorption require knowledge of photon trajectories at finite distances from the BH. As a particular example, one may consider an accretion disk around a black hole, a case which appears generically in astrophysics for BHs ranging in size from stellar-mass cases appearing as X-ray binaries up to Supermassive BHs observed as Active Galactic Nuclei. The intense thermal heating resulting from the collisional process in these cases can yield copious high-energy emission, while the densities are potentially large enough in some regions that both emission and absorption must be considered to properly describe the full radiation transfer present. In this work, we consider the full trajectory of photons that traverse the equatorial plane of a Kerr BH starting and ending at spatial infinity. We note that there is an immediate generalization for photons both emitted and absorbed at finite radius, $r$, as such photons superimpose onto trajectories that start and end at spatial infinity, and thus may be characterized by the parameters defined in figure 1.

Calculations of light deflection in general have also been an active area of research, frequently within the context of understanding gravitational lensing and its potential observational signatures. The complete solution for null geodesics in stationary BH spacetimes are expressed in terms of elliptic integrals (see [14, 15] for expressions involving Schwarzschild spacetimes, and [16] for those involving Kerr). Orbits in Kerr metrics and variants, including Kerr-de Sitter spacetimes with a cosmological constant, Kerr-Newman spacetimes with an electric charge, or both, have been computed exactly by Kraniotis [17, 18, 19] in terms of generalized multivariable hypergeometric functions and Weierstraß elliptic functions. The numerical evaluation of these integrals and/or complicated mathematical functions places a limit on the computational efficiency of ray-tracing codes that require large numbers of trajectories per volume. Approximate schemes have been developed for both Schwarzschild [20] and Kerr geodesics [21, 22], but the tradeoffs between accuracy and computational efficiency are apparent. Here, we report an approximate approach that provides a highly accurate closed-form expression that is also computationally efficient.

Our paper is organized as follows: In Section 2, we review the basic equations governing null geodesics in Schwarzschild and Kerr spacetimes and establish notation, such as the critical impact parameter, $b_{c}$, and radius of closest approach, $r_{0}$, shown in figure 1. In Section 3, we discuss the current understanding of light bending in terms of the strong-field ($b\to b_{c}$ in figure 1) and weak-field ($b\to\infty$ in figure 1) limits as a function of the BH spin. Here we also review the necessary information from Paper 1 regarding the bending angle, provide new coefficients for the weak-field asymptotic expansion, introduce the closest-approach ($r\to r_{0}$ in figure 1) and far-distance ($r\to\infty$ in figure 1) limits in the radial direction, and provide a recursive formulation for the full asymptotic expansion in both of these limits. In Section 4, we use the method of asymptotic approximants to construct full analytic trajectories of photons that bridge the closest approach and far distance limits for given impact parameters, and in doing so, naturally bridge the strong and weak field limits. In Section 5, the asymptotic approximant is compared with the full numerical solution for the trajectory of light. We conclude the main text of the paper in Section 6, with a summary of key findings, possible future refinements, and a short discussion contrasting asymptotic approximants and asymptotic matching.  A provides an important symmetry relationship involving the photon trajectories that enables efficient calculations.  B provides useful series relations and the steps for developing recurrence relations for the series coefficients in both the asymptotic limits examined as well as the approximants. For convenience, C provides a summary of coefficients needed to determine the overall bending angle ($\alpha$ in figure 1) used in the construction of the approximant – these coefficients are largely taken from Paper 1 with a few noted refinements.

We describe our photon trajectory in polar coordinates with radius $r$ and azimuthal angle $\phi$ (see figure 1), defining the inverse radius $u\equiv r^{-1}$ for convenience. Note that $r$ (and hence $u$) is taken here to be dimensionless, scaled by the distance unit $GM_{\rm BH}/c^{2}$ where $G$ is the gravitational constant, $M_{\rm BH}$ is the BH’s mass, and $c$ is the speed of light. The BH spin parameter $a$ is dimensionless as well, defined by

where $J_{\rm BH}$ is the BH’s angular momentum. The dimensionless spin parameter has a magnitude that naturally varies from $a=0$, a Schwarzschild BH without any angular momentum, up to $|a|=1$, an extremal Kerr BH. Measuring BH spins observationally is difficult, but several groups have provided estimates. For Sgr $A^{*}$, the supermassive black hole in the center of the Milky Way galaxy, estimates of the spin parameter range from $|a|=0.52$ from periodicities seen in infrared flares [23] up to spins of $a\approx 0.995$ as derived from periodicities in X-ray flares [24]. Gravitational wave detectors like LIGO (the Laser Interferometer Gravitational-wave Observatory) and VIRGO will someday be able to constrain BH spins through their effects on the waveforms of merging binaries, but to date are only able to measure moderate misalignment of spins between the two BHs prior to merger, particularly in the source GW170104 [25], without tight constraints on individual BH spins themselves. In this work, we simply note that the range of physically motivated parameters extends from $|a|=0$ up to $|a|\sim 0.999$, with the upper limit still a matter of theoretical uncertainty. The case $|a|=1$, which we also consider here, is included primarily for mathematical rather than astrophysically motivated reasons, as an asymptotic limit rather than a case likely to be observed in nature.

The Kerr geodesic equation for equatorial plane trajectories (thus having polar angle of $\pi/2$) is given by [5, 8]

where $\phi_{0}$ is the angle of closest approach and $\phi$ limits to $-\pi$ at $u=0$, as per figure 1. The non-dimensional impact parameter, $b$, is defined in terms of the $Z$-component of the photon’s angular momentum, $L$, and its energy, $E$, as

The relationship between $\phi$, $r$, and $b$ is shown in figure 1. As indicated, $b$ represents the minimum (perpendicular) distance between the BH center and the photon’s unperturbed trajectory a photon would travel in the absence of gravitational curvature effects.

We assume all photons travel with time, $t$, in a counterclockwise direction in the equatorial plane, with increasing azimuthal angle, $d\phi/dt>0$ (as indicated by the direction of the arrows along each trajectory in figure 1). For cases where the photon trajectory is prograde, we define the spin to be in the $+Z$ direction (out of the page in figure 1), with $0<a\leq 1$; for retrograde, the spin is in the $-Z$ direction and $-1\leq a<0$. By integrating 1, we find the azimuthal angle takes the form

The denominator of the integrand in (2) has zeros that distinguish its solution; these are, in fact, integrable singularities except along the ICO. Both the quadratic term and the cubic polynomial $h$ (defined in (1)) have real-valued zeros, but there is one positive root of the cubic that occurs at a smaller value of $u$; it thus dominates the behavior of the integral. The only exception to this situation is shown in Paper 1 for the extremal case $a=1$, where the zeros in the quadratic become coincident with those of $h$ as the ICO is approached, leading to a singular asymptotic behavior in that limit (see Paper 1: Appendix). With these preliminary comments, we examine the nature of the cubic in more detail.

The cubic polynomial in $u$ that appears as $h(u;a,b)$ in (1) has three real roots, two of which are always positive and one negative. We denote the smaller positive root as $u_{0}$; this is the root referred to above that dominates the integral behavior. This quantity, defined as $u_{0}=1/r_{0}$, is the the largest value of $u$ that a photon beginning from large distance ($u\rightarrow 0$) can achieve before reaching the closest approach along its trajectory (occurring at $r_{0}$ in figure 1), at which point it will begin to recede from the black hole. We note that if we convert $h(u;a,b)$ to $h(1/r;a,b)$ in (1) and multiply through by $r^{3}$, we arrive at a cubic in $r$. The distance of closest approach, $r_{0}$, may be found in terms of $a$ and $b$ by solving the cubic equation $h(1/r;a,b)$=0:

Eq. (3) implies that $r_{0}$ and $b$ are monotonically related. Thus, for a given spin, photon trajectories can be parameterized by either $r_{0}$ or $b$.

In Paper 1, we considered the relationship between a given impact parameter, $b$, and the total bending angle, $\alpha$, both situated at infinite radial distances from the BH as shown in figure 1; it was possible to do so without the details of the trajectory, according to the formula

Note that the integral in (4) evaluates to $\pi/2$ for a straight line trajectory in the absence of a gravitating source, yielding $\alpha=0$, consistent with the geometry in figure 1.

As written, (1) and all further equations involving $\phi$ are valid in a portion of the domain in figure 1, corresponding to $\phi\in(-\pi,\phi_{0})$, for which $d\phi/dr<0$ (i.e., $d\phi/du>0$) and $r\in(\infty,r_{0})$ (i.e., $u\in(0,u_{0})$). As shown in Paper 1, the remaining region corresponding to $\phi^{\prime}\in(\phi_{0},\alpha)$, and for which $d\phi^{\prime}/dr>0$ (i.e., $d\phi^{\prime}/du<0$), is described by (1) (replacing $\phi$ with $\phi^{\prime}$) and with a right-hand side of $-f$. This observation reveals that the full trajectory is symmetric about $r=r_{0}$ (or $u=u_{0}$). To map out the trajectory of a photon, the procedure taken is to prescribe an impact parameter, $b$ (and thus, $r_{0}$), at $r=\infty$, which corresponds to $u=0$. At any finite distance $r$ , the local angle $\phi$ may be determined from (2) where $u\leq u_{0}$. Once $\phi$ is obtained, the local ($X,Y$) coordinates of the trajectory are extracted as $X=r\cos(\phi)$ and $Y=r\sin(\phi)$ in the region $\phi\in(-\pi,\phi_{0})$. The remaining portion of the trajectory in $\phi^{\prime}\in(\phi_{0},\alpha)$ can be mapped out using the symmetry condition of (23) in A.

In this work, we construct approximate solutions for the trajectory of photons as they pass a Kerr BH, as described by (1), such that the procedure described above can be carried out in an accurate but efficient manner. Note that although a photon trajectory is characterized by the impact parameter that is defined at spatial infinity, photons themselves do not necessarily need to start or end at spatial infinity along a given trajectory. The parameter space for such a trajectory in $(r,\theta)$ is encompassed within 4 limits, shown in figure 2a and each described in Section 3. Here we introduce the required notation and a convenient alternative coordinate system. In order to normalize the domain in $u$, we define a new quantity

to use as an integration variable, chosen so that the upper bound of the integral (2) satisfies $y=1$ in the limiting case and $0\leq y\leq 1$ in general. We also adopt the same convention as [7], and normalize the impact parameter (shown in figure 2b) as

where $b_{c}$ (shown in figure 2a) corresponds to the limiting case of infinite bending angle, given by

which follows from finding the $b$ value that minimizes $r_{0}$ in (3)²²2This occurs at the value of $b$ at which the argument of the arccosine in (3) is $-1$, and is in general a solution of a cubic equation.. It also follows from (3) that the critical minimum separation, $r_{c}$ (shown in figure 2a), for a given spin, $a$, is

where one may verify that $r_{0}=r_{c}$ when $b=b_{c}$ in (3). The effect of (5) and (6) is to map the infinite physical domain of figure 2a to the unit square of figure 2b. It is convenient to utilize the latter coordinate system and $b^{\prime}$ notation when ultimately displaying results. However, for purposes of notational clarity in our development, we will leave integrals and their expansions in terms of $b$.

Applying the variable transformation (5) to the integral (2), the azimuthal angle takes the following form in terms of $y$:

which is linked to the bending angle, $\alpha$, considered in Paper 1 (and described by (4)) via the expression

where $\alpha=\alpha(a,b)$ and $\phi_{0}$ (the angle of closest approach) are both shown in figure 1. Our goal here is to determine approximate formulae for $\phi$ in terms of the three parameters $a$, $b$, and $y$. As in Paper 1, we will treat the BH spin $a$ solely as a parameter. In contrast to that work, in which the impact parameter $b$ was generally treated as a variable, here we use $y$ as the primary expansion variable for our approximants. There is physical motivation for this choice: a trajectory with fixed $b$ but variable $y$ represents the evolution of a photon through the BH space-time, as shown in the schematic in figure 1 (noting that $y=r_{0}/r$).

We now proceed to determine limiting asymptotic expressions for the azimuthal angle, $\phi$ (shown in figures 1 and 2), in various limits. Note that expressions are derived in the mapped domain (figure 2b) in what follows, but the reader is may wish to refer back to the physical domain (figure 2a) to obtain a clear picture of the physical limits being examined. Also note that in what follows, the BH spin, $a$, is treated solely as a parameter, and is assumed fixed.

In Paper 1, we restricted attention to the overall deflection of the photon, $\alpha$, for a given scaled impact parameter, $b^{\prime}$ (related to the usual impact parameter according to (6), as shown in figure 1); this is related to the value of the azimuthal angle along the boundary $y=1$ (i.e., $\phi_{0}$ according to (10)) as shown in figure 2b. Our approximant for the bending angle $\alpha$ was formed by combining results from the strong-field limit, corresponding to $b^{\prime}\rightarrow 0$, and the weak-field limit, corresponding to $b^{\prime}\rightarrow 1$, and using asymptotic approximants to bridge the gap for intermediate values of $b^{\prime}$. As shown in figure 2b, there are two additional limits that correspond to the left and right-hand side of the domain, where now $b^{\prime}$ is treated as a parameter. In particular, the far-distance limit, for which $y\rightarrow 0$, corresponds to cases where the BH trajectory begins at spatial infinity but ends arbitrarily far from the BH as well. The opposite case is the closest-approach limit, for which $y\rightarrow 1$ and the photon trajectory ends arbitrarily close to the pericenter passage. Note that for large values of $b$ (i.e., $b^{\prime}\to 1$), the closest-approach may lie very far from the BH, while for trajectories with small values of $y$, the strong-field limit may apply for photons that are never actually close to the BH. Regardless of the choice of the BH spin $a$, the azimuthal angle remains finite over the entire 2-dimensional phase space except for the strong-field, closest-approach case $b^{\prime}=0$ and $y=1$ (upper-right corner of figure 2b), which corresponds to the critical trajectory for which a photon asymptotically approaches the innermost circular orbit and never escapes the BH.

While expansions can be developed for any of the 4 limits described above, the behavior of the azimuthal angle in these various limits differs in important ways. An effective canvassing of the space is needed to construct an asymptotic approximant, and is accomplished using the particular expansions given in the following subsections. Since the expansion variable, $y$, is the same as the integration variable in (9) for the far-distance and closest-approach limits, we choose to traverse the space of figure 2b by slicing horizontally. This is also physically motivated, since each horizontal line in figure 2b corresponds to a photon trajectory, as shown in figure 2a. The zeroth-order term of the closest-approach expansion is $\phi=\phi_{0}(a,b^{\prime})$ (also shown in figure 2), which is not known exactly but an accurate expression (asymptotic approximant developed in Paper 1) is given in Section 3.1 below, which includes new terms (higher-order in $b^{\prime}$) not disclosed previously. All higher-order terms of the closest-approach expansion are determined exactly in Section 3.2. The zeroth-order term of the far-distance expansion is $\phi=-\pi$, as shown in figure 2 (a or b). All higher-order terms of the far-distance expansion are determined exactly in Section 3.3. The approximant of Section 4 is formed using the expansions given in Sections 3.1 through 3.3, and thus naturally incorporates all 4 limits of figure 2b.

To overcome the issue of divergent series such as those given in (18) and (17), there are a host of convergence acceleration (or “re-summation”) techniques available that rely only on the original series itself (e.g. Padé approximants, continued fractions, Euler summation, etc.) [27]. While such methods typically lead to an implementation improvement compared with the original series, global accuracy is not always guaranteed and the “best” re-summation technique is not always obvious [28, 29, 30]. In particular, any finite representation of an infinite power series can be analytically continued in a number of ways. Additionally, summation techniques, such as Padés (i.e. rational functions), can be effective because they approximate the pole singularities responsible for series divergence. Of course, if the singularity responsible for the divergence is not a pole, the technique is not often as effective. And, since the function we are trying to approximate is often unknown (only having its divergent series expansion to work with), it can thus be difficult to choose an efficient re-summation method a priori.

Asymptotic approximants constrain the analytic continuation of a series derived in one limit so that it can approach an asymptotic behavior in a different limit; an overview with examples is given in [13]. The accuracy of the approximant may be improved by including additional terms in the power series representation used in one of the limits. The approach is validated by forming a convergent sequence of approximants as additional terms are added. In the current work, the far-distance and closest-approach asymptotic expansions, given respectively by (18) and (17) are used to construct an approximant (see figure 2).

The current problem mirrors a problem in thermodynamics that has been solved successfully using asymptotic approximants [11, 12], namely bridging a low-density power series at one end with a non-integer power law at the other end (the thermodynamic critical point). The non-integer asymptotic behavior renders the rational-function form of standard Padés innefective. In that problem, only knowledge of the leading-order non-integer behavior is required to construct a uniformly accurate asymptotic approximant. Similarly, in the current problem, we have a regular power series in the far-distance limit and a non-integer asymptotic expansion in the closest-approach limit. And akin to the thermodynamics problem, we find that only the leading-order behavior of the closest-approach expansion is required to form a uniformly accurate approximant for any spin $a$ and most values of $b$ (higher-orders needed as $b\to b_{c}$), as is demonstrated in what follows (in Section 5).

An asymptotic approximant for the trajectory is constructed by satisfying the asymptotic expansions for the far-distance and closest-approach limits for a given impact parameter $b$ shown in figure 2a, corresponding to $b^{\prime}$ in figure 2b. The strong-field and weak-field limits are implicitly incorporated since their asymptotic forms are consistent with the closest-approach and far-distance limit expressions for small and large values of the impact parameter. An approximant that satisfies the far-distance limit (18) to $N^{\mathrm{th}}$-order and the $K$-term non-integer expansion in the closest-approach limit (17) is given by

where $\phi_{0}$=$(\alpha-\pi)/2$, $\tilde{\tilde{g}}_{0}=\tilde{g}_{0}-\phi_{0}$, $\tilde{\tilde{g}}_{n>0}=\tilde{g}_{n>0}$ (defined in (18)), and $C_{n}$ coefficients are given in (17). For all approximant curves generated here, $\alpha$ (imbedded in $\phi_{0}=(\alpha-\pi)/2$) is computed using the bending angle approximant of Paper 1 taken to 5^th order, provided in Section 3.1 and with additional details in C. The steps for obtaining the expressions for $A_{n}$ and $T_{n}$ in (19) are given in B.3. The approximant (19) matches the closest-approach expansion (as $y\to 1$) asymptotically to $(1/2+K)^{th}$-order (choosing $K\geq 0$) and has an expansion about $y=0$ that is exactly the far distance series to $N^{\mathrm{th}}$-order. Note that one may also set $K=-1$ to remove the $K$ series from (19). This form of the approximant enforces only the zeroth-order closest approach limit $\phi_{0}$, but still matches the functional form of higher-order terms. For $K=-1$, the approximant (19) reduces to

The reduced approximant given by (20) may be used to compute accurate photon trajectories within most of physical parameter space (i.e., that of figure 2b). This is demonstrated in the cases presented in Section 5, and guidance is then provided for when one might use the higher-order approximant given by (19) for $K\geq 0$.

The approximant (19) is compared with the numerical solution of (9) in figures 6 through 10. All trajectories in the figures are constructed using either the approximant (solid curves) or the numerical solution ($\circ$’s), with the aid of the symmetry relations given in A. Figure 6 shows the trajectory of photons around an extremal ($a=1$) BH, presented in two coordinate systems: the physical $X$-$Y$ plane (see figure 1), used in figures 6a and 6b; and the mapped coordinate system (see figure 2b), used in figures 6c and 6d. In what follows we present results in both coordinate systems, as the former is physically appealing, while the latter is more sensitive to the accuracy of the approximant.

In figure 6a trajectories for $b^{\prime}$= 0.1, 0.5, and 0.7 are produced by the reduced $K=-1$ approximant given by (20) using $N$=6. Note that, even for the $b^{\prime}$=0.1 curve where photons pass thrice around the BH, the $K=-1$ approximant is indistinguishable from the numerical solution on the scale of the figure. This is almost the case for the even closer trajectory ($b^{\prime}$=0.05) shown in figure 6b, where photons pass around the black hole six times. Here, the $K=-1$ approximnant (20) (dotted line) becomes misaligned with the numerical solution (see figure inset) as the photons move radially closer to the BH. This issue is corrected by including more information from the closest-approach limit, as shown by approximant (19) using $K=1$ (the solid curve of figure 6b, see inset). While the difference in approximants is difficult to see on the scale of figure 6b, the mapped coordinate system (showing the same information) in figure 6c clearly shows the error in the $K=-1$ approximant, visible between $y\approx 0.7$ and $y\approx 0.95$. Here, the $K=-1$ approximant (20) (labeled A6) uses only information from $\phi_{0}$ (indicated on figure) in the closest-approach limit and the far-distance series (18) using $N$=6 terms (labeled F6 on figure). The higher-order approximant (19) (also shown in figure 6c, labeled A15) uses information from the closest-approach expansion (17) (labeled C1) taken to 3/2-order (i.e. $K=1$) and the far-distance series (18) using $N$=15 terms (labeled F15 on figure); this curve is indistinguishable from the numerical solution on the scale of the figure. This approximant is extended over a range of $b^{\prime}$ to produce the surface shown in figure 6d.

Note that the value of $N$ is different for each of the two approximants discussed above (and used in figure 6). This is attributed to the fact that the sequences of approximants we have constructed are in fact divergent past a certain order, and there is an optimal number of terms to obtain a best approximation in each case. Although both the closest-approach and far-distance expansions are known to infinite order, they are divergent series, each having a radius of convergence (dependent on $a$ and $b^{\prime}$) within the physical domain $y\in[0,1]$ as shown by the plots in figure 3 (Section 3.2) and figure 4 (Section 3.3). That said, the singularities responsible for the radius of convergence are non-physical, as they are not resident in the integrand of (9) in the interval $y\in[0,1]$. The effect of this is that the approximant has an ‘optimal truncation’ [27], $N=N_{\mathrm{opt}}$, such that the $(N_{\mathrm{opt}}+1)^{\mathrm{th}}$ term leads to a smaller contribution to the approximant than both lower and higher-order terms (i.e. the remainder series diverges). This feature is illustrated in figure 7 for an extremal ($a=1$) BH with $b^{\prime}=0.1$ using the reduced $K=-1$ approximant (20); here, $N_{\mathrm{opt}}=6$, as indicated in figure 7b.

As additional terms are included from the closest-approach series (increasing $K$) in approximant (19), $N_{\mathrm{opt}}$ changes. The effect of $K$ on optimal truncation and achievable accuracy in the approximant is shown in figure 8, where the relative error is plotted for $K=-1,~0,~1$, and 2 for an extremal ($a=1$) BH with $b^{\prime}=0.1$. In each subfigure (representing different $K$), $N$ is taken up to its optimal truncation, as discussed above. Note that for each $K$, as $N$ increases, the error decreases in some region of $y$ (prior to optimal truncation) but the norm of the error never goes below $O$(10^-4). In Section 6, we suggest possible ways in which this persistent error may be reduced.

In order to demonstrate its versatility, the approximant is shown for $a=0.95$ and $a=0$ in figures 9 and 10, respectively. In each figure, 4 plots are provided: (a) the photon trajectory, using the reduced $K=-1$ approximant (20); (b) the $\phi$ surface in the $y$ vs. $b^{\prime}$ plane (using, if required, the higher order approximant (19), i.e., $K\neq 1$); (c) a cross section of the $\phi$ surface at $b^{\prime}=0.5$, showing a comparison between the $N$-term far-distance (F$N$) and $K$-term closest-approach (C$K$, or $\phi_{0}$ for $K=-1$) expansions that are used as an input to the approximant, where $N$ is taken to be $\leq N_{\mathrm{opt}}$ such that the approximant is accurate on the scale of the figure; and (d) the relative error of approximant (20) for increasing $N$ (up to optimal truncation $N_{\mathrm{opt}}$) at a fixed $K$. The values of $a$ and $b^{\prime}$ for these figures are chosen to show the range of the approximant’s usage while preserving accuracy. Note that accuracy is maintained for other values of $a$ and $b^{\prime}$ provided that additional terms are used from the closest-approach expansion as $b^{\prime}\to 0$ and additional terms are used from the weak-field limit (in $\phi_{0}$) as $a\to 1$, as discussed in Paper 1.

We now provide some guidance on using approximant (19) and its reduced form (20). For plotting trajectories of any spin $a$ in the range $b^{\prime}\in[0.1,1]$, such as those shown in figures 6a, 6b, 9a, and 10a, one may use the reduced $K=-1$ approximant (20) for simplicity’s sake. If used for visualization purposes, these trajectories will be indistinguishable from the numerical solution, with an error of (at most) $O$($10^{-3}$). For closer trajectories ($b^{\prime}<0.1$), or in any case where one wishes to reduce local error, we recommend using the more general approximant (19), which includes corrections from the closest-approach expansion. One may then determine the optimal truncation by following the procedure illustrated in figure 7b, which does not require knowledge of the exact solution.

In regards to computational expense, the authors have found the closed-form solution given by (19) to be, in general, an order of magnitude faster than the numerical evaluation³³3The numerical evaluation of (9) is done using the “quad” command in MATLAB, which uses an adaptive recursive Simpson’s rule [31] of (9), when generating the same number of data points. Furthermore, the authors found that time-savings from using less terms (either in $K$ or $N$) in the approximant are negligible. We recommend using as many terms as are required for the desired accuracy, as it will be faster than solving the exact elliptic integrals numerically.

Analytical expressions that describe the full trajectory of photons propagating in the equatorial plane of a Kerr black hole are obtained using asymptotic approximants. These have potential use in future ray-tracing efforts and radiation transport numerical projects, particularly those studying accretions disks around black holes and their observable electromagnetic emission. The expressions obtained provide accurate trajectory predictions for arbitrary spin and impact parameters, and provide significant time advantages compared with numerical evaluation of the elliptic integrals that describe photon trajectories. The asymptotic approximants provided here are accurate closed-form expressions that bridge the weak-field (large impact parameter), strong-field (near-critical impact parameter), far-distance (large radius), and closest approach (smallest radial distance from the black hole) limits. To that end, asymptotic expansions are derived for the azimuthal angle in the far-distance and closest-approach limits, and new coefficients are reported for the bending angle in the weak-field limit.

While an optimal truncation of the approximants provided here can be determined to maintain accuracy, future work may focus on incorporating the non-physical singularities responsible for the divergence of the far-distance and closest approach series. One idea would be to use a modified Padé approach to incorporate these singularities as poles, while maintaining the correct asymptotic behavior in all regions listed above. This would involve a matrix inversion in the process of computing the approximant coefficients, whereas now their computation is solely recursive. That said, there are efficient methods for computing this particular inversion [32] and this would still likely require less computational time than the full numerical solution of the elliptic integrals.

We close this paper by distinguishing the approach taken here, i.e., asymptotic approximants, from that of asymptotic matching used in singular perturbation theory. In the later approach, a region of overlap must occur in which two asymptotic expansions are valid, and a systematic method is used to obtain an expression that is uniformly valid over the whole domain. The key is that, in asymptotic matching, the physical domain is completely encompassed and described by overlapping asymptotic expansions. Although asymptotic approximants do connect two asymptotic expansions in different regions of a domain, it is not necessary that the original asymptotic expansions overlap, and there can be a gap in which both expansions are not valid (see, for example, figure 5). By forming a convergent sequence (via analytic continuation) of approximants, the gap region is well approximated where both asymptotic series can fail.