Above is a figure of the projected emission rings for different (solid), along with lines of constant (dashed). The double-valued nature of the parameterization is here visualized.
A degenerate parameterisation of the image plane may be constructed, since every photon observed emitted from the accretion disc may be associated with a in the bundle of photons coming from , and that there is an implicit map .
Using a change of variable substitution, the observed flux is
where the partial differential term is a Jacobian.
Cunningham transfer functions were originally defined in Cunningham's paper, however the precise quoting in subsequent papers by other authors sometimes drop a normalizing term. Cunningham justifies the original definition
This form of has been chosen so that its numerical value is nearly independent of and .
This is presumably to solve numerical issues.
The transfer function methods have been popularized by a number of works in computational astrophysics, including Dauser et al. (2010), in which a Green's function formulation espoused additional convenience. For that reason, I use the definition presented in Dauser et al. (2010),
Note that , like , is double-valued except for extremal redshifts coming from . Integrating over the appropriate limits
presents a number of computational problems, both in determining , and in avoiding singular values as where the Jacobian diverges.
Throughout this section, we will assume a geometrically thin disc in the equatorial plane. The methods are trivially generalizable to other geometries by projecting all radial components of the disc into the equatorial plane.
Let us begin by discussing methods for calculating the transfer functions. The first step is calculating the parameterization. An algorithm may be as follows:
For a target emission radius , assume the circular ring on the disc may be projected as a star-shaped boundary on the image plane, given by pairs of . Since this region is star shaped, one may express
and find for a given for which the geodesic intersects the disc at radius . The set of give which may be used to determine the set of corresponding to the chosen emission radius, from which , are found, and used to construct the parameter.
The computational difficulty of this algorithm is in accurately determining , and sampling enough to have a good estimate of , .
This problem is simple to solve when rephrased as a root finding problem:
function find_ξ_for_rₑ(metric, u, disc, rₑ, ϑ) # define a map ξ → (rₑ - r) on the disc 𝔉 = ξ -> begin α = ξ * cos(ϑ) β = ξ * sin(ϑ) point = integrate_single_geodesic(metric, u, disc, α, β) r = point.u2 rₑ - r end # solve with Order0 : https://juliamath.github.io/Roots.jl/dev/reference/#Roots.Order0 ξ_sol = Roots.find_zero(𝔉, 2rₑ + 20.0) return ξ_sol end
Note that the measure to find the root is and not, for example, the other way around. This is done as photons hitting the event horizon have increasing , despite their projection into the disc is decreasing. Subtracting might therefore present two unique roots for small , which is undesirable, or at the very least prevent bracketing limits from having opposite signs.
The choice of initial value here is arbitrary to scale with emission radius. A more informed prior will improve the performance of the algorithm.
Accurately estimating the extrema of may be achieved similarly, using the resulting to determine , and then pass this function to a non-linear solver to determine which map to extremal . The only problem here is that each step of the non-linear solver must root-find, which quickly becomes a costly problem.
Other authors (Bambi et al. (2017) and Abdikamalov et al. (2020)) use a coarse version of this to effectively binary search for the minima and maxima close to . Our method is to use a cubic spline interpolation over the calculated in the previous step, and then root-find over the derivative of the interpolation. This is both fast and surprisingly accurate in Julia, due to automatic-differentiation-enabled interpolations libraries like PumasAI/DataInterpolations.jl, and automatic-differentiation backends like JuliaDiff/ForwardDiff.jl.
We apriori know that the minima and maxima of will be close to , however the domain of puts the minima and maxima at the edges of the domain, which may be difficult to optimize. We shift the domain of to , such that the minima and maxima are close to and respectively, in the middle of our domain.
With as few as 10 knots, the cubic interpolation already achieves a good estimate of the extrema. In practice, we use a minimum of around 20 knots. These are currently always equally spaced, however different sampling methods may be explored in the future. With the current setup, the error relative to taking the extrema of 2000 equally spaced points scales in the following manner:
Also shown (brown) are the errors of taking the extrema of knots without interpolation, as a baseline from which the interpolation improves.
The periodicity exhibited is due to aliasing effects when subdividing the interpolation knots equally. We can expect the error to be for our chosen 20 knots relative to taking the extrema of 2000 points.
The code used to find the maximum and minimum is
∂(f) = x -> ForwardDiff.derivative(f, x) function interpolate_extremal(y, x, x0) interp = DataInterpolations.CubicSpline(y, x) x̄ = Roots.find_zero(∂(interp), x0) x̄, interp(x̄) end
which makes use of an initial guess
x0, i.e. or .
Bambi et al. (2017) describe a method for determining the Jacobian calculations involved in , and were kind enough to send me their code to examine when mine was failing. Their method makes use of
to calculate the differential terms by offsetting and , and numerically calculating the gradient at those points using effectively a 1 order finite differencing method.
Any finite differencing stencil may be applied to some central , and this is something we have exploited with help of JuliaDiff/FiniteDifferences.jl:
function jacobian_∂αβ∂rg(metric, u, disc, α, β; diff_order = 5) # map impact parameters to r, g 𝔉 = ((α_, β_),) -> begin point = tracegeodesic(metric, u, disc, α_, β_) g = redshift(metric, point) # return r and g [point.u2, g] end cfdm = FiniteDifferences.central_fdm(diff_order, 1) J = FiniteDifferences.jacobian(cfdm, 𝔉, ([α, β])) |> first abs(inv(det(J))) end
We calculate the Jacobian with respect to and not purely for optimization reasons: in our code it is cheaper (from a memory + compute perspective) to defer calculating and until after all of the and had been calculated for a given , as then and could be discarded (see Gradus source code). The Jacobian is then rescaled with
All of the required components for calculating the transfer functions , as defined in eq. (6), are now known. Below are a number of illustrative for for different viewer inclinations, as labelled in the curves, for two different black hole spins.
The magnitude of the transfer functions changes dramatically for different viewing angles and spins at , since both of these alter the projected velocity of the accretion disc: changing the viewing angle changes the component of the velocity parallel to the line of sight, whereas changing the spin both alters the radius of the ISCO and the photon momentum. In the above figures, this is enough to reorder the ISCO transfer function contributions of different viewing angles.
Integrating the transfer functions to produce flux profiles as a function of energy requires a degree of attention. Starting from eq. (8), Dauser et al. (2010) introduce their Green's function formalism with the substitution
equating the source emission to a delta function at a specific energy . For convenience, use that
for the given specific energy , and then substitute into the flux integrand and evaluate the delta by integrating over . I had to remind myself that delta functions have the property
and so we gain a factor in the integrand, but without loss of generality one can set . We also perform a single variable substitution of to keep the resulting expression familiar, and pick up corresponding factors:
Since we want to find the flux at different energies, we split the integral into discrete bins of , equivalent to :
Continuing with this integral notation betrays the nature of the calculation we are attempting, so lets start digging into the details and walk through them step-by-step. Principally, there are three problems to address
Handling the double-valued nature of by splitting the integral into two branches.
Faithfully interpolating the transfer functions over and .
Avoiding singular (divergent) values of at when integrating.
The first point requires separating the transfer function in continuous domains of by splitting the transfer functions at the points where and . There are a number of possible ways to do this, but our approach is to group the knots into an upper and lower branch by walking through all of the points sequentially, and then interpolate both branches over . This has the drawback that points close to are required for the interpolation to be faithful, but knots seems to be sufficient enough in practice. In the context of the integral, we split and integrate as normal.
The second and third point may be addressed together.
Integrating the transfer functions according to eq. (16) is non-trivial. Bambi et al. (2017) and Dauser et al. (2010) describe integrating with respect to first, and then over . The temptation is then to interpolate over , and then marginalize the integral. However we, like Dauser et al. (2010), approach this the other way round. The motivation for this comes purely from numerical stability – the integration algorithms seemed to converge faster when approached this way.
The breadth and depth of numerical integration is beyond this post, but our method makes use of a interval based integration method (trapezoidal integration), and an adaptive Gauss-Kronrod quadrature scheme, as these methods harmonize well with the nature of the problem. Gauss-Kronrod additionally has the benefit that some interval , the integrand is never evaluated directly at or , allowing us to avoid the singularities in at extremal . The requirement that the integrand be smooth is also accounted for in our cubic spline interpolation.
Dauser et al. (2010) handle the divergent points by examining the limits of the transfer functions
concluding that the total integrand of eq. (16) diverges as , and that by assuming over the small integration interval, the analytic solution goes as
They qualify this by , and use this approximation in the regions of and , where the normalizing factors are determined from and .
In practice, the diverging region is finite and also related to unstable Jacobian values, and consequently , even using Gauss-Kronrod integration, does not entirely vanish. Dauser et al. (2010) use for their approximation, using Romberg adaptive trapezoid integration for . Our method can set , and effectively ignore the divergent region entirely, however we chose to keep the Dauser et al. (2010) approximation over this small until we have had the opportunity to explore the errors in more depth.
Trapezoidal integration evaluates the integrand at discrete locations and uses a weighted sum to calculate the integral. This permits lazy evaluation, where the weight may be calculated prior to the integrand, and the same weight used for any integrand over the same interval.
Let us illustrate this with a code snippet from the integration:
function integrate_drdg✶(ε, transfer_functions, radii, g_grid; N=10) # pre-allocate output flux = zeros(Float64, length(g_grid)) minrₑ, maxrₑ = extrema(radii) # create N knots of interpolation interp = interpolate_over_radii(transfer_functions, N) # build fine radial grid for trapezoidal integration fine_rₑ_grid = build_rₑ_grid(minrₑ, maxrₑ, 1000) |> collect for (i, rₑ) in enumerate(fine_rₑ_grid) # wrap integrand integrand = wrap_interpolations(interp) # trapezoidal integration weight if i == 1 Δrₑ = (fine_rₑ_grid[i+1]) - (rₑ) elseif i == lastindex(fine_rₑ_grid) Δrₑ = (rₑ) - (fine_rₑ_grid[i-1]) else Δrₑ = (fine_rₑ_grid[i+1]) - (fine_rₑ_grid[i-1]) end # all radial factors weight = Δrₑ * rₑ * ε(rₑ) # integrate each bin in the redshift grid for j in eachindex( (g_grid[1:end-1])) glo = g_grid[j] ghi = g_grid[j+1] flux[j] += integrate_bin(integrand, rₑ, glo, ghi) * weight end end return flux end
A grid of is created for the trapezoidal integration, in this case with 1000 intervals. The weight for each interval may be calculated once, and then used for each transfer function over between
ghi, i.e. the limits of the energy domain we wish to calculate the flux for.
All of the integrand terms applying to are scoped in a closure during
wrap_interpolations. The precise implementation is unimportant, and can be seen in the Gradus.jl code repository as mentioned.
interpolate_over_radii function requires closer examination: so far, we have a discrete set of transfer functions separated into upper and lower branches, but now we wish to evaluate the transfer functions for arbitrary , and therefore we interpolate over along each of the different values of , i.e. the knots.
integrate_bin uses the interpolations to evaluate the integrand. Keeping with Dauser et al (2010), we handle the divergent region separately:
function integrate_bin(integrand, rₑ, lo, hi; h = 2e-8) gmin = integrand.gmin(rₑ) gmax = integrand.gmax(rₑ) # ensure we don't go out of bounds glo = clamp(lo, gmin, gmax) ghi = clamp(hi, gmin, gmax) intensity = 0.0 # no bin width, i.e. out of limits for this transfer function if glo == ghi return intensity end g✶lo = g_to_g✶(glo, gmin, gmax) g✶hi = g_to_g✶(ghi, gmin, gmax) if (g✶lo < h) || (g✶hi > 1-h) # ... handle edge integration ... end res, _ = QuadGK.quadgk(integrand, glo, ghi) intensity += res return intensity end
The bulk of the work is performed here by JuliaMath/QuadGK.jl. We will discuss the edge approximation in just a moment. For the full source code, see the Gradus.jl source.
For posterity, I want to note here that my initial implementation had the transfer function interpolations mapping , and integrating over . I spent many weeks debugging this implementation, before realizing that I had implicitly performed a change of variable in the integrand, and was missing a factor . The current implementation of the code uses interpolations that map to mitigate this.
The normalizing factor for the "edge" region of the integration is a little more complex in practice than just using scaling values of as described in Dauser et al. (2010). Indeed, their Fortran code does something a little more akin to:
function integrate_edge(integrand, h, lim, gmin, gmax, low::Bool) if low gh = g✶_to_g(h, gmin, gmax) a = √gh - √lim else gh = g✶_to_g(1 - h, gmin, gmax) a = √lim - √gh end 2 * integrand(gh) * √h * a end
The factor 2 comes from integrating the diverging terms that go as , whereas I know is related to the bin width, but am unsure of its motivation. Dauser et al. (2010) have an additional factor , which I assume is a Jacobian term, but including any of these factors actually seems to make my normalization worse when is artificially increased. This function is non-critical, rarely invoked, and mostly returns values extremely close to zero when it is – consequently I have very little incentive to investigate the accuracy of the normalizing terms at the moment, though I certainly will revisit this at some point!
Emerge from the swamps of detail, let us look at the fruits of our labour: we compare the line profiles obtained using Gradus.jl to Dauser et al. (2010)
I took the liberty of bootstrapping
relline (and all of the XSPEC model library) for Julia with the incredible work of JuliaPackaging/BinaryBuilder.jl. All of this is available in the SpectralFitting.jl package I am working on.
using SpectralFitting d = GeometricThinDisc(0.0, 400.0, π / 2) u = [0.0, 1000.0, deg2rad(40), 0.0] m = BoyerLindquistAD(M=1.0, a=0.998) redshift = ConstPointFunctions.redshift # maximal integration radius maxrₑ = 50.0 # emissivity function ε(r) = r^(-3) # g grid to do flux integration over gs = range(0.0, 1.2, 500) _, flux = lineprofile(gs, ε, m, u, d, redshift_pf = redshift, maxrₑ = maxrₑ) # transform to observed energy energy = gs .* 6.4 # invoke relline flux_relline = invokemodel( energy, XS_Relline( outer_r = FrozenFitParam(maxrₑ), θ_obs = FrozenFitParam(rad2deg(u)), limb = FrozenFitParam(0), ), )
Producing the comparison plot:
This is quite good but there are obvious noise contributions that need to reduced (discussed in the next section).
On my university 2021 M1 Mac running with 4 threads the whole line profile calculation takes approximately 16 seconds, with almost all of that time spent on calculating the transfer functions. Part of the intention with Gradus.jl is being able to generate novel spectral models on the fly, and with this performance, the entire
relline table model data could be generated at the same resolution in approximately half an hour on my laptop. If we would therefore be interested in, for example, different spacetimes, accretion disc geometries, or more, and have priors on the parameter ranges we are interested in, smaller, targeted tables could be generated within the time it takes to make a cup of coffee!
Speaking of other spacetimes, our implementation is generic – merely passing in a different metric and redshift point function allows us to produce line profiles for tests of relativity:
m = JohannsenPsaltisAD(M=1.0, a=0.6, ϵ3=2.0) redshift = interpolate_redshift(interpolate_plunging_velocities(m), u) # ...
Here we use the Johannsen-Psaltis metric from Johannsen & Psaltis (2011) for a variety of deformation parameters with otherwise the same disc and observer configuration as before:
My initial implementation for calculating the transfer functions works well for , but after this deviates from the Dauser et al. (2010) calculated table, and indeed even begins to produce increasingly large transfer function magnitudes for larger radii, where the colloquial lower branch miraculously becomes the upper branch. My understanding is that this issue related to how the redshift is calculated at large impact parameters , , corresponding to a wide field-of-view.
I had previously noticed that assuming for the photon energy at the observer adds a vignetting error to the redshift, since an observer at distance feels only a rough approximation of flat spacetime. Instead calculating properly at the observer gives energies slighty different from , approximately 1 part in 100 if I remember correctly scaling with . This difference would further be exacerbated by the root finding, mapping, and Jacobian calculator.
This should then be resolved in two ways: either use the proper observer energy calculations, or simply move the observer back, scaling with , keeping , small. I've opted for the latter, to keep the redshift equations consistent with those in Cunningham and Fabian et al. (1997). I will investigate the source of this deviation further at a later date.
Here, the left panel does not move the observer back as a function of , whereas the right panel does.
Another caveat is in the finite difference stencil: depending on the sensitivity of the transfer function Jacobian at high or , these methods sometimes require different finite differencing orders, with no single value that works well for all cases. Using a different algorithm (central vs forward, etc.) could help here, but I am tempted instead to use automatic-differentiation, since the SciML/DifferentialEquations.jl solvers support dual number types. This is something that will likely change in my implementation in the future, however for now we just dynamically change the stencil order.
For example, using a finite stencil of order is too coarse at large emission radii, where the differences in redshifts over an emission ring can be very small with at , versus at for . This is reflected in numerical noise in the calculated transfer functions:
Left panel is central finite difference stencil with order , right panel is order .
Since most emissivity models for assume some dependence, the terms at large radii are suppressed, and the noise in the functions is less important. I believe it is still worth mentioning, to illustrate how fiddly these calculations can be.
|||C T Cunningham (1975). The effects of redshifts and focusing on the spectrum of an accretion disk around a Kerr black hole. APJ, 202:788-802. NASA ADS|
|||Misner, C. W., Thorne, K. S., Wheeler, J. A., and Kaiser, D. I., Gravitation. 2018. NASA ADS|
|||Dauser, T., J. Garcia, J. Wilms, M. Böck, L. W. Brenneman, M. Falanga, K. Fukumura, and C. S. Reynolds. “Irradiation of an Accretion Disc by a Jet: General Properties and Implications for Spin Measurements of Black Holes.” Monthly Notices of the Royal Astronomical Society 430, no. 3 (April 11, 2013): 1694–1708. DOI.|
|||Bambi, Cosimo, Alejandro Cárdenas-Avendaño, Thomas Dauser, Javier A. García, and Sourabh Nampalliwar. “Testing the Kerr Black Hole Hypothesis Using X-Ray Reflection Spectroscopy.” The Astrophysical Journal 842, no. 2 (June 15, 2017): 76. DOI.|
|||Abdikamalov, Askar B., Dimitry Ayzenberg, Cosimo Bambi, Thomas Dauser, Javier A. Garcia, Sourabh Nampalliwar, Ashutosh Tripathi, and Menglei Zhou. “Testing the Kerr Black Hole Hypothesis Using X-Ray Reflection Spectroscopy and a Thin Disk Model with Finite Thickness.” The Astrophysical Journal 899, no. 1 (August 14, 2020): 80. DOI.|
|||Johannsen, Tim, and Dimitrios Psaltis. “A Metric for Rapidly Spinning Black Holes Suitable for Strong-Field Tests of the No-Hair Theorem.” Physical Review D 83, no. 12 (June 7, 2011): 124015. https://doi.org/10.1103/PhysRevD.83.124015.|