ABSTRACT (original article): Hotspots, often characterized as pointlike emissions, frequently appear near black holes with significantly enhanced luminosity compared to the surrounding accretion flow. Notably, such hotspots are regularly observed near the black hole at the center of the Milky Way. Light rays emitted from these sources follow complex trajectories around the black hole before reaching distinct locations on the observer’s image plane. Precisely resolving both direct emissions and their higher-order images—despite the latter’s intensity suppression—is essential for extracting detailed spacetime information, including the black hole’s mass, spin, and inclination angle. To improve the accuracy and efficiency of hotspot modeling, we develop a forward ray tracing method based on the analytic integral solution of Kerr geodesics, leveraging conserved quantities. Our approach traces geodesics from a given emission point near the black hole to a distant observer, effectively capturing multiple images with a tailored parametrization scheme for root-finding. By perturbing these geodesics, we map finite-size emissions to distinct regions on the image plane, enabling the quantification of image shapes and amplification rates. This method not only enhances the identification of strongly lensed photons from black holes but also enables efficient spacetime tomography and hotspot localization, leveraging observations from the Event Horizon Telescope and its upcoming next-generation upgrades. CITATION (original article): Lihang Zhou, Zhen Zhong, Yifan Chen, Vitor Cardoso (2025), Forward ray tracing and hot spots in Kerr spacetime, Phys. Rev. D 111, 064075. https://doi.org/10.1103/PhysRevD.111.064075
Github: https://github.com/AuroraDysis/KerrP2P
Github: https://github.com/AuroraDysis/KerrP2P
KerrP2P is a software designed for forward ray tracing in Kerr spacetime. It is specifically tailored to efficiently calculate multiple null geodesics between designated “source” and “observer” points, locate apparent positions of the corresponding images, and quantify their shapes. Detailed information can be found in the paper Forward Ray Tracing and Hot Spots in Kerr Spacetime by Lihang Zhou, Zhen Zhong, Yifan Chen, and Vitor Cardoso.
This software contains a Python/C++ part and a Mathematica part. The Python/C++ part is aimed at high-precision contour plot calculation to locate the multiple images. The Mathematica part (i.e., this notebook) takes the conserved quantities generated by the former as input, visualizes the geodesics in 3D, and calculates the positions and shapes of the resulting images.
Def: Ray
Def: Ray
In this section define the Ray function, which calculates the final direction of a geodesic based on its starting location and initial direction encoded by conserved quantities.
Prerequisites
Prerequisites
Ray
Ray
Def: FindImage
Def: FindImage
In this section we define the FindImage function, which solves the imaging conditions to give the conserved quantities of a geodesic that reaches the observer, and then calculates multiple quantities that describe the trajectory of the geodesic, including the half-orbit number , polar turning number , azimuthal winding angle , travel time , apparent position (α, β) of the resulting image and so on.
This function provides two kinds of input guess values:
1. together with an string opt=“in” or “out” specifying the sign of (“in” for and “out” for )
2. (λ, q)
Due to the fact that higher-order images exponentially approach the critical curve, it would require more precise guess values of for larger image levels, and therefore in that case the (rc, |d|) input is preferred.
(=and≡(mod2π))
θ
f
θ
o
ϕ
f
ϕ
o
n
m
ϕ
f
t
f
This function provides two kinds of input guess values:
1.
(rc,|d|)
log
10
d
d<0
d>0
2. (λ, q)
Due to the fact that higher-order images exponentially approach the critical curve, it would require more precise guess values of
(λ,q)
log
10
In[]:=
(*findimagebyinputting(rc,d)asinitialguess*)FindImage::opt="opt should be in or out";FindImage::νrs="νrs should be -1 or +1";FindImage::νθs="νθs should be -1 or +1";Options[FindImage]={"WorkingPrecision"->30};FindImage[a0_,rs0_,θs0_,ro0_,νrs_Integer,νθs_Integer,θo0_,ϕo0_,{rc1_?NumericQ,lgAbsd1_?NumericQ},opt_String,OptionsPattern[]]:=Module{sol,rcsol,dsol,m,nhalf,λ,η,r1,r2,r3,r4,w0,root$f,opt$sign,α,β,err,ray,a,rs,θs,ro,rc0,lgAbsd,θo,ϕo,θf,ϕf},If[opt=!="in"&&opt=!="out",Message[FindImage::opt];Return[FindImage::opt]];If[νrs=!=+1&&νrs=!=-1,Message[FindImage::νrs];Return[FindImage::νrs]];If[νθs=!=+1&&νθs=!=-1,Message[FindImage::νθs];Return[FindImage::νθs]];opt$sign=If[opt==="out",1,-1];{a,rs,θs,ro,θo,ϕo,rc0,lgAbsd}=SetPrecision[{a0,rs0,θs0,ro0,θo0,ϕo0,rc1,lgAbsd1},OptionValue[WorkingPrecision]];root$f[rc_?NumericQ,lgd_?NumericQ]:=Withrayf=Rayrcda,rs,θs,ro,νrs,νθs,rc,opt$sign,WorkingPrecision->OptionValue[WorkingPrecision]+10,CalculateFinalTime->False,rayf["θf"]-θo,Sin;err=Catch[sol=FindRoot[root$f[rc,lgd],{{rc,rc0},{lgd,lgAbsd}},WorkingPrecision->OptionValue[WorkingPrecision]];];If[err=!=Null,Return[err]];{rcsol,dsol}=rc,opt$sign/.sol;ray=Rayrcd[a,rs,θs,ro,νrs,νθs,rcsol,dsol,WorkingPrecision->OptionValue[WorkingPrecision]+10];{θf,ϕf,m,nhalf}={ray["θf"],ray["ϕf"],ray["m"],ray["nhalf"]};{λ,η}=rcd$to$λη[a,rcsol,dsol];w0=Cos[θpm[a,λ,η][[2]]];{α,β}=conserve$to$ob[a,λ,η,θf,νθs];{r1,r2,r3,r4}=radial$potential$roots[a,λ,η];<|"a"->a,"λ"->λ,"η"->η,"rs"->rs,"θs"->θs,"ro"->ro,"τo"->ray["τo"],"θf"->ray["θf"],"ϕf"->ray["ϕf"],"tf"->ray["tf"],"θo-θf"->θo-ray["θf"],"Sin"->Sin,"m"->m,"nhalf"->nhalf,"Mod[ϕf, 2π]"->Mod[ϕf,2π],"α"->α,"β"->β,"νrs"->νrs,"νθs"->νθs,"rc"->rcsol,"d"->dsol,"r1"->r1,"r2"->r2,"r3"->r3,"r4"->r4,"w0"->Cos[θpm[a,λ,η][[2]]],"radial_integrals"->ray["radial_integrals"],"angular_integrals"->ray["angular_integrals"]|>;
lgd
10
rayf["ϕf"]-ϕo
2
lgd
10
m
(-1)
ϕf-ϕo
2
ϕf-ϕo
2
In[]:=
(*findimagebyinputting(λ,q)asinitialguess*)FindImageλq::νrs="νrs should be -1 or +1";FindImageλq::νθs="νθs should be -1 or +1";Options[FindImageλq]={WorkingPrecision->30};FindImageλq[a0_,rs0_,θs0_,ro0_,νrs_Integer,νθs_Integer,θo0_,ϕo0_,{λ1_?NumericQ,q1_?NumericQ},OptionsPattern[]]:=Module{sol,λsol,qsol,m,nhalf,rc,d,r1,r2,r3,r4,w0,root$f,α,β,err,ray,a,rs,θs,ro,λ0,q0,θo,ϕo,θf,ϕf},If[νrs=!=+1&&νrs=!=-1,Message[FindImage::νrs];Return[FindImage::νrs]];If[νθs=!=+1&&νθs=!=-1,Message[FindImage::νθs];Return[FindImage::νθs]];{a,rs,θs,ro,θo,ϕo,λ0,q0}=SetPrecision[{a0,rs0,θs0,ro0,θo0,ϕo0,λ1,q1},OptionValue[WorkingPrecision]];root$f[λ_?NumericQ,q_?NumericQ]:=Withrayf=Rayληa,rs,θs,ro,νrs,νθs,λ,,WorkingPrecision->OptionValue[WorkingPrecision]+10,CalculateFinalTime->False,rayf["θf"]-θo,Sin;err=Catch[sol=FindRoot[root$f[λ,q],{{λ,λ0},{q,q0}},WorkingPrecision->OptionValue[WorkingPrecision]];];If[err=!=Null,Return[err]];{λsol,qsol}={λ,q}/.sol;ray=(*Rayrcd*)Rayληa,rs,θs,ro,νrs,νθs,λsol,,WorkingPrecision->OptionValue[WorkingPrecision]+10;{θf,ϕf,m,nhalf}={ray["θf"],ray["ϕf"],ray["m"],ray["nhalf"]};{rc,d}=λq$to$rcd[a,λsol,qsol];w0=Cos[θpm[a,λsol,][[2]]];{α,β}=conserve$to$ob[a,λsol,,θf,νθs];{r1,r2,r3,r4}=radial$potential$roots[a,λsol,];<|"a"->a,"λ"λsol,"η"->,"q"->qsol,"rs"->rs,"θs"->θs,"ro"->ro,"τo"->ray["τo"],"θf"->ray["θf"],"ϕf"->ray["ϕf"],"tf"->ray["tf"],"θo-θf"->θo-ray["θf"],"Sin"->Sin,"m"->m,"nhalf"->nhalf,"Mod[ϕf, 2π]"->Mod[ϕf,2π],"α"->α,"β"->β,"νrs"->νrs,"νθs"->νθs,"rc"rc,"d"d,"r1"->r1,"r2"->r2,"r3"->r3,"r4"->r4,"w0"->Cos[θpm[a,λsol,][[2]]],"radial_integrals"->ray["radial_integrals"],"angular_integrals"->ray["angular_integrals"]|>;
2
q
rayf["ϕf"]-ϕo
2
2
qsol
2
qsol
2
qsol
m
(-1)
2
qsol
2
qsol
ϕf-ϕo
2
ϕf-ϕo
2
2
qsol
Def: 3D Geodesic
Def: 3D Geodesic
In this section we first define a numeric code that integrates the geodesic equation to obtain the parameterized trajectory of a geodesic. Then we define a function used to demonstrate the geodesic in 3D.
Numeric geodesic
Numeric geodesic
3D Demo
3D Demo
Def: Shape
Def: Shape
This section is for calculating the shape of the observed image, based on the perturbative method used in the original paper.
Evaluate the ℳ matrix
Evaluate the ℳ matrix
Here we use both semi-analytic and numeric methods to evaluate the ℳ matrix characterizing the map between the source position and the image plane.
The semi-analytic expressions are derived by working out the relevant partial derivatives directly from the integral formalism, including properly treating the integrand and possible turning points. The results are expressed as definite integrals that will be evaluated numerically, so this method is dubbed as “semi-analytic”. The process of deriving it is quite lengthy and will not be shown here.
On the other hand, the numeric method is much more straightforward——numeric derivative by introducing a small perturbation ϵ to the parameters. We have checked that the two methods lead to consistent results.
The semi-analytic expressions are derived by working out the relevant partial derivatives directly from the integral formalism, including properly treating the integrand and possible turning points. The results are expressed as definite integrals that will be evaluated numerically, so this method is dubbed as “semi-analytic”. The process of deriving it is quite lengthy and will not be shown here.
On the other hand, the numeric method is much more straightforward——numeric derivative by introducing a small perturbation ϵ to the parameters. We have checked that the two methods lead to consistent results.
Semi-analytic derivatives
Semi-analytic derivatives
SA ℳ
SA ℳ
Numeric ℳ
Numeric ℳ
Consistency check
Consistency check
Display the shape
Display the shape
Examples
Examples
In this section we present examples to show the multiple functions of this Mathematica code.
Due to the high nonlinearity of the geodesic equation, it’s recommended to first solve the imaging conditions using contour plots (as described in the original paper using the python code in https://github.com/AuroraDysis/KerrP2P ), and then plug the values of conserved quantities into this Mathematica code to investigate geodesic visualization and image-plane position and shape.
(=and≡(mod2π))
θ
f
θ
o
ϕ
f
ϕ
o
Ex1: N=1 image
Ex1: image
N=1
This example appeared in the “Basic Usages” part of the Jupyter notebook tutorial.
In[]:=
a=;rs=10;θs=;θo=2.5;ϕo=-5;ro=1000;νrs=-1;νθs=-1;(*Settheguessedvaluesof(λ,q)tofindaimageandextractdetailedinformation.*)λ$seed=-2.9;q$seed=5.3;data=FindImageλq[a,rs,θs,ro,νrs,νθs,θo,ϕo,{λ$seed,q$seed},"WorkingPrecision"->80];If[StringQ@data,Print[data],Grid[Transpose@{Keys@data,Values@data},Frame->{All,All}]]
8
10
89.9π
180
3D demonstration of the geodesic . Increase smax to extend the geodesic .
Image plane
Image shape
Ex2: Higher-order image in the paper
Ex2: Higher-order image in the paper
3 D demonstration of the geodesic . Increase smax to extend the geodesic .
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Forward ray tracing and hot spots in Kerr spacetime
by Lihang Zhou
Wolfram Community, STAFF PICKS, April 16, 2025
https://community.wolfram.com/groups/-/m/t/3444948
by Lihang Zhou
Wolfram Community, STAFF PICKS, April 16, 2025
https://community.wolfram.com/groups/-/m/t/3444948