© The Authors 2026

1. Introduction

The Event Horizon Telescope (EHT) has produced images of the supermassive black holes Messier 87* (hereafter M87*) and Sagittarius A* (hereafter Sgr A*), providing an unprecedented laboratory for testing accretion models and studying the effects of light propagation in general relativity (EHTC 2019b, 2022b, 2024a). Additionally, observations by the GRAVITY collaboration have constrained the mass of Sgr A* to M = 4.3 × 10⁶ M_⊙  ± 0.25% (GRAVITY Collaboration 2022), establishing a crucial reference for event-horizon-scale tests of space-time, such as black hole spin measurements and potential deviations from general relativity.

A major challenge when performing precision space-time tests from black hole images is disentangling the plasma effects of accretion from the pure general relativity signatures of the underlying space-time. The combination of finite observational resolution, strong interstellar scattering along the Galactic Center line of sight, and the short dynamical timescale of Sgr A* (approximately 30 minutes at the innermost stable circular orbit for a Schwarzschild black hole) makes this separation particularly difficult. Consequently, any attempt to estimate the dimensionless spin parameter, a_* = Jc/GM², must rely on identifying a signature in one of the less plasma-dependent observables.

Among the most distinctive features of black hole space-times is their ability to trap photons in unstable yet bound orbits (Bardeen et al. 1972). In the case of a Kerr black hole, the shape and size of these orbits depend solely on the black hole’s mass and spin. The region where these bound orbits reside is known as the photon sphere. As photons escape these unstable orbits and reach a distant observer, they appear as a narrow “photon ring” on the image plane; they can asymptotically approach the so-called critical curve depending on the number of windings the photon undergoes around the black hole before arrival. These features have been extensively studied in recent works, both theoretically (Claudel et al. 2001; Johannsen & Psaltis 2010; Gralla et al. 2019; Gralla & Lupsasca 2020; Kocherlakota et al. 2024) and in connection with observations (Johnson et al. 2020; Gralla et al. 2020; Wong 2021; Zhou et al. 2025; Zhang et al. 2025). The Black Hole Explorer (BHEX) mission (Johnson et al. 2024; Lupsasca et al. 2024; Galison et al. 2024) has set the precise measurement of the photon ring as a key science goal.

In this work we focused on a specific effect of black hole lensing: the secondary images of localized sources, or “hot spots”, a subcategory of photon rings. Unlike a full ring, the secondary image remains compact and forms a crescent. This provides a unique opportunity to distinguish direct from lensed photons in the image, enabling the application of the mathematical framework developed for null geodesics. We studied two key observational effects: the time lag between primary and secondary photons (n = 0 and n = 1) and the difference in position angle on the screen (ΔPA), with a particular focus on the latter. Specifically, we derived a new analytic approximation for calculating δPA (where δPA denotes the analytic estimate, while ΔPA refers to the observational measurement), propose a method for disentangling the inclination angle (i) from ΔPA measurements, and establish an empirical relation connecting ΔPA to black hole spin. Definitions of the different terms used in this paper for ΔPA are provided in Table 1.

This method is entirely achromatic, as it depends solely on the geometric properties of the system rather than specific emission mechanisms. While the appearance of a single hot spot may differ between, for example, radio and near-infrared (NIR) wavelengths, the geometric relations underlying our method remain valid in both regimes.

Observations of Sgr A* support the presence of emission from localized regions, as the source exhibits regular flaring events across the electromagnetic spectrum on characteristic timescales of ∼20 hours. During these events, NIR radiation can increase by 1–2 orders of magnitude (GRAVITY Collaboration 2020a), while the brightest X-ray flares reach peaks of ∼600 times the quiescent state (Haggard et al. 2019). Additionally, GRAVITY has detected astrometric motion of the Sgr A* brightness centroid and polarimetric signatures associated with flares (GRAVITY Collaboration 2018, 2023). At radio frequencies, observations with the Atacama Large Millimeter/submillimeter Array (ALMA) have been used to link similar physical effects to corresponding signatures in the light curves of linear polarization (Wielgus et al. 2022a,b; EHTC 2022c). Proposed next-generation improvements to the EHT have specifically targeted improved dynamical imaging sensitivity for the Galactic Center due to the presence of these rapidly evolving events (Doeleman et al. 2023).

Numerous studies have modeled this localized emission (e.g., Broderick & Loeb 2006; Trippe et al. 2007; Hamaus et al. 2009; Tiede et al. 2020; Gelles et al. 2021; Vos et al. 2022; Vincent et al. 2024), and some have sought to connect simulations to observations for model classification (GRAVITY Collaboration 2020b,c; Ball et al. 2021; Aimar et al. 2023; Yfantis et al. 2024a,b; Antonopoulou & Nathanail 2024). As these models increasingly reproduce observational data, a stronger theoretical foundation becomes essential.

The broader concept of hot spots has been linked to physical processes in magnetically arrested disks (MADs; Narayan et al. 2003), which exhibit flux eruption events resembling the flaring activity of Sgr A* (Dexter et al. 2020; Scepi et al. 2022). These events can generate flux tubes – evacuated regions in the accretion flow where electrons become highly accelerated and heated via magnetic reconnection while traveling in strong vertical magnetic fields (Porth et al. 2021; Ripperda et al. 2022). Simulations by Najafi-Ziyazi et al. (2024) suggest that such structures produce signatures consistent with observations. As MADs have sub-Keplerian orbital profiles (Begelman et al. 2022; Conroy et al. 2023), the hot spots from flux eruptions should also appear sub-Keplerian. We expect such phenomena to be most prominent in the submillimeter regime, as was shown in Yfantis et al. (2024b), where a sub-Keplerian velocity for the hot spot was strongly favored from the fitting of ALMA data. In this case we would also expect the bulk of the emission to come from the equatorial plane, since the eruptions are stronger there. Our analysis focuses more on this scenario because equatorial orbits are simpler, but also because all the near-term future experiments capable of imaging in the required resolution are in the radio regime.

Another possible explanation is the formation of plasmoids via magnetic reconnection in magnetized disks (e.g., Ripperda et al. 2020; Aimar et al. 2023; Vos et al. 2024). In this case the modeling could get more complicated, with off-equatorial orbits (spiraling inward or outward) and modular velocities, as plasmoids can detach from the bulk of the disk much more easily. These effects cannot be expressed in this analysis and could significantly impact the observables. Notably though, plasmoids are expected mostly in the NIR and higher frequencies, with velocities above the Keplerian limit, differentiating them from the simpler case. In the future we plan to extend the method to include all this rich phenomenology.

An alternative phenomenon that could appear as bright “spots” are the so-called shock waves (appearing as spiral arms), as shown in Conroy et al. (2023). They are known to reside in the equatorial plane, but it is not clear whether they are strong enough that their secondary imprint would survive.

Lastly, it should be possible for two bright spots to appear simultaneously, although it is not clear at what cadence. Regardless, we could use various gauges to differentiate this from the case of a secondary image, such as the time delay between the appearance of the two images, the flux difference, and their orbital properties.

Given the growing scientific progress in both observational and modeling efforts, we anticipate that a direct image or time-resolved movie of a secondary image will become a reality in the near future. In this paper we present a straightforward and computationally efficient approach to interpreting such observations – contrasting with the full analytic treatment of photon rings – with the goal of constraining black hole spin. To this end, we have constructed a library of over 900 models with varying key system parameters: the spot’s orbital radius (r_hs), the orbital velocity, the black hole spin (a_*), and the inclination angle (i).

In Sect. 2 we summarize the key mathematical results from previous studies. In Sect. 3 we present our numerical scheme and the model library we used to study the problem. In Sect. 4 we introduce our intermediate findings: a simple and highly accurate (sub-degree) approximation for δPA_0 → 1, (as opposed to the n → ∞ that is typically considered), along with a simplified formula that has a growing deviation with spin from the sub-degree level up to 5° for a_* = 0.99. We then show that the mean value of ΔPA over a full orbit of the hot spot is insensitive to i, similar to recent results in Walia et al. (2025). This result favors simple analyses of average deflections, removing the notion of a time-dependent orbital phase and the observer’s inclination and position angle in downstream fitting. Lastly, we introduce an empirical relation derived from our model library and evaluate its accuracy in reproducing simulation results. We then apply this relation to three sets of mock observations, incorporating reasonable uncertainties for ΔPA, orbital period (P), and r_hs. We show that in the case of observations that align with our simulation assumptions, our method can constrain the black hole spin to ±0.3 at the 2σ confidence level, but this fitting should be mainly viewed as a proof of concept; a dedicated fitting analysis should be performed in the future. In Sect. 5 we conclude by discussing the implications and potential limitations of this approach for future observations.

2. Mathematical prior

The goal of this section is to provide a short summary of the mathematical background governing the effects of lensing around black holes. We provide this information as a means for comparison and as a natural path toward an empirical relation that stems from simulations and which incorporates the equations presented here.

We used Boyer-Lindquist coordinates (t, r, θ, ϕ) and natural units, i.e., G = c = M = 1, so that both time and distance are measured in [M]. The line element describing the space-time is

$$\begin{array}{cc}& d{s}^2=-(1-\frac{2r}{\mathrm{\Sigma }})d{t}^2-\frac{2a_{\ast }r{sin}^2\theta }{\mathrm{\Sigma }}dtd\varphi +\frac{\mathrm{\Sigma }}{\mathrm{\Delta }}d{r}^2+\mathrm{\Sigma }d{\theta }^2+\hfill \\ \hfill & (r^2+a_{\ast }^2+\frac{2a_{\ast }^{2}r{sin}^2\theta }{\mathrm{\Sigma }}){sin}^2\theta d{\varphi }^2,\hfill \end{array}$$(1)

where a_* is the dimensionless angular momentum per unit mass and

$$\begin{array}{c}\hfill \mathrm{\Sigma }(r,\theta )=r^2+a_{\ast }^2{cos}^2(\theta ),\mathrm{\Delta }(r)=r^2-2r+a_{\ast }^2.\end{array}$$(2)

We considered a stationary hot spot and show solutions for two quantities: the time lapse (δt) between two subsequent images (n, n+1) and the difference in their position angle on the sky (δPA) for an observer at infinity, on the symmetry axis. The solution holds only for n → ∞, implying that the parameters are evaluated at r_crit. In this paper we omit the complete solution, which contains various steps and different cases; it can be found in numerous papers (e.g., Beckwith & Done 2005; Gralla & Lupsasca 2020; Wong 2021). Instead we provide the solutions (similarly to Gralla & Lupsasca 2020) as

$$\begin{array}{cc}\hfill \delta t& =\frac{2}{\sqrt{{\stackrel{\sim }{b}}^2-a_{\ast }^2}}[{\stackrel{\sim }{r}}_0^2\left(\frac{{\stackrel{\sim }{r}}_0+3}{{\stackrel{\sim }{r}}_0-1}\right)K\left(\frac{a_{\ast }^2}{a_{\ast }^2-{\stackrel{\sim }{b}}^2}\right)\hfill \\ \hfill & -2a_{\ast }^2{E}^{\prime }\left(\frac{a_{\ast }^2}{a_{\ast }^2-{\stackrel{\sim }{b}}^2}\right)]\hfill \end{array}$$(3)

and

$$\begin{array}{c}\hfill \delta \text{PA}=\pi +\frac{2a_{\ast }}{\sqrt{{\stackrel{\sim }{b}}^2-a_{\ast }^2}}\left(\frac{{\stackrel{\sim }{r}}_0+1}{{\stackrel{\sim }{r}}_0-1}\right)K\left(\frac{a_{\ast }^2}{a_{\ast }^2-{\stackrel{\sim }{b}}^2}\right),\end{array}$$(4)

where

$$\begin{array}{c}\hfill {\stackrel{\sim }{r}}_0=1+2\sqrt{1-\frac{a_{\ast }^2}{3}}cos[\frac{1}{3}arccos\left(\frac{1-a_{\ast }^2}{{(1-\frac{a_{\ast }^2}{3})}^{3/2}}\right)]\end{array}$$(5)

is the critical radius, denoting the bound photon orbit for the face-on observer (i = 0°). The tilde on the rest of the parameters indicates that they are calculated on the critical radius. Then $\stackrel{\sim }{b}$, the apparent critical curve on an observer’s screen in units of M, is given by

$$\begin{array}{c}\hfill \stackrel{\sim }{b}=\sqrt{\frac{{\stackrel{\sim }{r}}_0^3}{a_{\ast }^2}[\frac{4\mathrm{\Delta }({\stackrel{\sim }{r}}_0)}{{({\stackrel{\sim }{r}}_0-1)}^2}-{\stackrel{\sim }{r}}_0]+a_{\ast }^2}.\end{array}$$(6)

Lastly, K and E′, in Eqs. (3) and (4), represent the first and second complete elliptical integrals respectively, which can be found in math libraries, such as scipy.

For the case of a moving hot spot there is an extra degree of complexity, since the spot changes position while the photons move to the secondary image, resulting in a necessary correction for δPA (Eq. (4)). In the simplest case of a circular orbit on the equator it is given by

$$\begin{array}{c}\hfill \delta {\text{PA}}_{\text{orb}}=\delta \text{PA}-\delta t\mathrm{\Omega },\end{array}$$(7)

where Ω is the angular velocity of the spot.

As stated, and further demonstrated in Appendix A, these equations describe the problem only for higher-order images (n → ∞). In this paper, however, we are interested in the secondary image (n = 1), so we need to quantify how much the solution deviates in this case. To our knowledge, there is no published solution for such a case; hence, we present a heuristic approximation in Sect. 4.

3. Numerical method

In this section we describe our numerical setup, built inside the general relativistic radiative transfer code ipole (Mościbrodzka & Gammie 2018). We simulated emission from a moving hot spot, closely resembling that implemented in Yfantis et al. (2024a,b). The hot spot lives in circular orbits on the equatorial plane (mostly expected in radio flares), and it has a Gaussian profile peaking at the center for its key parameters: number density (n_e), electron temperature (Θ_e) and magnetic field strength (B_field). Its size is dictated by the number density that has a 1σ radius of 1.5M and then a steep cutoff, creating a crisp spot size. While the spot model is defined in Boyer–Lindquist coordinates, the ray-tracing of photons is done in Kerr-Schild coordinates. The resolution is set at 256 × 256 pixels.

Our library consists of models with different values for a variety of parameters (r_hs,  Ω_coef,  i, and a_*), where r_hs is the distance of the spot from the center of the black hole and the Ω_coef parameter controls the orbital velocity as

$$\begin{array}{c}\hfill P=\frac{2\pi }{{\mathrm{\Omega }}_{\mathrm{coef}}}(r_{\text{hs}}^{1.5}+a_{\ast }).\end{array}$$(8)

The Ω_coef parameter provides a straightforward way to simulate sub-Keplerian orbits, highly expected in MADs, or even slower orbits associated with spiral shock waves. Appendix A provides proof that Ω_coef maps correctly to the velocity in the image domain. The parameter grid is summarized in Table 2.

During the simulation we recorded values of the centroids of the two images, and we masked all contributions to light rays from n > 1. A simulation that includes the n = 2 image, showcasing that its contribution is completely negligible, is presented in Appendix A. Given the centroids, we defined ΔPA as

$$\begin{array}{c}\hfill \mathrm{\Delta }PA=\text{atan2}(Y_{{\text{cen}}_2}{X}_{\text{cen}}-X_{{\text{cen}}_2}{Y}_{\text{cen}},X_{\text{cen}}{X}_{{\text{cen}}_2}+Y_{\text{cen}}{Y}_{{\text{cen}}_2}),\end{array}$$(9)

where X_cen, Y_cen and X_cen2, Y_cen2 refer to the centroid of total intensity of primary and secondary images, respectively. It is worth mentioning that the correct ΔPA of a moving source can only be estimated using the slow light approach and that is what we did throughout the study. Using the fast light approximation, the second term of Eq. (7) cannot be manifested, so it is effectively a stationary hot spot, leading to large deviations¹.

We also measured the arrival time of all photons to the observer, which allowed us to calculate the maximum time difference (Δt) between the primary and secondary image as the difference between the first and last photon to reach. In Fig. 1 we present two example snapshots from the simulation domain with different parameters. The imaging happens after the secondary image has appeared, meaning 5 − 15 minutes after the hot spot has emerged. The middle and right panels display the evolution of ΔPA and Δt over a full orbit. ΔPA is calculated as the angle starting at n = 0 and ending at n = 1, counterclockwise (see Fig. 1), while Δt as the difference in time between the last photon from n = 0 and the first from n = 1. Both quantities exhibit significant variation with increasing inclination, highlighting the need for a method for disentangling inclination effects from other influences. This is presented in Sect. 4.

Fig. 1. Left: Snapshot from a hot spot movie, with resolution 256 × 256 (as in the entire library), where ΔPA is calculated from the primary to the secondary counterclockwise. Middle: Curve of ΔPA measurement with respect to time for a full period. Right: Time lag (difference in arrival times) between the first and second image. The gray vertical line in the right panels denotes the snapshot’s time in the left panels.

4. Results

This section includes two intermediate results necessary for the construction of an empirical relation: (1) an analytic approximation for δPA_0 → 1 and (2) the independence of the mean ΔPA (from a full period) with respect to inclination, similarly to the independence of the median value shown in Walia et al. (2025). Table 1 summarizes the notation of the approximations discussed here. Next we present our empirical relation that covers shortcomings from the analytic relations in order to be used as a means of inference, and lastly we show three examples of its application using mock observations of ΔPA.

4.1. Insightful analytic approximations for δPA_n → 1

To address this problem, we created three models with spin parameters a_* = 0,  0.5,  0.95, a smaller spot size (spot radius of 1 M), and a higher resolution (1024 × 1024 pixels). The spot is stationary and r_hs is set to the critical radius for all models.

First we tested the simulation results against the δPA solution. While both methods give ΔPA = 180° for the non-spinning case, for a_* = 0.5,  0.95 the analytic solution gives δPA = 216° ,263° while the simulations produce ΔPA = 198° ,218°, respectively. The reason is that the trajectories are sufficiently different for n = 0 and n = 1 that the critical parameters do not apply exactly. Thus, the complete integral (K) in Eq. (4) should be calculated as incomplete instead (K_inc). These are connected via

$$\begin{array}{c}\hfill K(x)=K_{\text{inc}}\left(\frac{\pi }{2}|x\right),\end{array}$$(10)

where the argument in K_inc is an angle in the first part and the calculated quantity in the second. So for π/2 the incomplete becomes complete.

Via trial and error, we found that the correct angle for an incomplete integral to account for the primary-secondary (0 → 1) image solution (δPA_0 → 1) is π/(π + a_* + 0.3) such that

$$\begin{array}{cc}& \delta {\text{PA}}_{0\to 1}=\pi +\frac{2M{a}_{\ast }}{\sqrt{{\stackrel{\sim }{b}}^2-M^2{a}_{\ast }^2}}\left(\frac{{\stackrel{\sim }{r}}_0+M}{{\stackrel{\sim }{r}}_0-M}\right)\times \hfill \\ \hfill & K_{\text{inc}}\left(\frac{\pi }{(\pi +a_{\ast }+0.3)}|\frac{a_{\ast }^2}{a_{\ast }^2-{\stackrel{\sim }{b}}^2}\right).\hfill \end{array}$$(11)

The reason behind the a_* dependence inside K_inc is that the deviation of the two formulas (δPA, δPA_0 → 1) increases for high spins, as the photons are traveling a larger angle (δPA increases with increasing spin), and that needs to be accounted for. Another effective and significantly simpler way to calculate the deviation in δPA_0 → 1 is by taking the expression for δPA (Eq. (4)) and removing the factor of 2 from the second term. We call this δPA_{f2,  0 → 1}, and it is given as

$$\begin{array}{c}\hfill \delta {\text{PA}}_{f2,0\to 1}=\pi +\frac{a_{\ast }}{\sqrt{{\stackrel{\sim }{b}}^2-a_{\ast }^2}}\left(\frac{{\stackrel{\sim }{r}}_0+1}{{\stackrel{\sim }{r}}_0-1}\right)K\left(\frac{a_{\ast }^2}{a_{\ast }^2-{\stackrel{\sim }{b}}^2}\right).\end{array}$$(12)

This works nearly exactly (sub-degree deviations) for spins up to ∼0.75 with deviations of just a couple degrees above that (5° for a_* = 0.99). A naive but perhaps useful intuition for this approximation is that the n = 0 photons do not feel the twist of the space-time in their way in, only as they move toward the back of the black hole, and so the deflection between n = 1 and n = 0 is only half of the effect in the universal regime.

Despite these simplifications, the analytic problem remains challenging due to two key generalizations: the motion of the hot spot and the treatment of a general emission radius r_hs, which is not fixed at the critical radius. We therefore proceeded with a broad numerical parameter study to identify a simpler empirical relation.

Fig. 2. Distribution (violins) of ΔPA values during a full period for different inclination values. The crosses denote the average values per period and the dashed lines the value at i = 0° for a fixed a*.

4.2. Inclination dependence

To study the variations in ΔPA with inclination, we created distributions (violins) of all ΔPA values, sampled every dt = 5M over a full period, for various models. These are shown in Fig. 2, where the models are characterized by two fixed parameter sets (r_hs = 4, 8 and Ω_coef = 0.5, 0.75), three spin values (a_* = −0.9, 0, 0.9), and five inclination angles in the range [0, 60]. The two parameter sets (r_hs,  Ω_coef) are shown in separate panels, while all spin values are overlaid within each panel. The vertical axis represents ΔPA, and the horizontal axis denotes the inclination angle.

The distributions shown in Fig. 2 exhibit significant variations in both ΔPA values and their shapes depending on the parameters. However, a clear pattern emerges: the mean ΔPA for each model, marked by crosses and connected with dashed lines, remains nearly constant across inclinations. This holds especially well for inclinations up to 30° and for larger values of r_hs, like in r_hs = 8, where all average ΔPAs are within one degree. Even in the case of r_hs = 4, the deviation at i = 60° from i = 0° is 5°. In the case of Sgr A* inclination is constrained by various observations to be i ≲ 30° (e.g., GRAVITY Collaboration 2018, 2020b; EHTC 2022d, 2024b; Wielgus et al. 2022b; Levis et al. 2024; Yfantis et al. 2024b,a), making this particularly effective. This is possible only in the case of a full orbit, since in smaller fractions of a period, specialized modeling would be required to understand the exact phase that is observed, for example comparing velocity and ΔPA curves of different inclinations. So, given a full-period dataset, the mean ΔPA can be reliably used for masking inclination effects. Furthermore, the spread of ΔPA values provides a potential means for estimating inclination itself. This important property enables the creation of a simple empirical relation.

4.3. Empirical relation for ΔPA

By analyzing a variety of the modeled ΔPA curves, and testing possible dependences using simple fitting python libraries, we found the following empirical relationship connecting three fundamental parameters of a hot spot around a black hole (r_hs,  Ω_coef,  a_*) with the observed angle difference between primary and secondary image (ΔPA):

$$\begin{array}{c}\hfill \mathrm{\Delta }PA{[}^{°}]={180}^{°}+a_{\ast }(40-r_{\text{hs}})-(\frac{1408-233a_{\ast }}{r_{\text{hs}}^2}+35){\mathrm{\Omega }}_{\mathrm{coef}},\end{array}$$(13)

where Ω_coef, can be expressed in terms of P,  r_hs and a_* by rearranging Eq. (8) (and using physical units) as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{coef}}=2\pi \frac{r_{\text{hs}}^{1.5}+a_{\ast }}{P[\text{M}]}=2.2248\frac{r_{\text{hs}}^{1.5}+a_{\ast }}{P[\text{min}]},\end{array}$$(14)

where the last expression includes mass of Sgr A* from GRAVITY Collaboration (2022), M = 4.3 × 10⁶ M_⊙.

There are three key improvements of Eq. (13) with respect to the deflection Eqs. (11) and (12):

1. The addition of the angular lag caused by the movement of the emitter. This is present in Eq. (7), but the complexity increases significantly with the inclusion of δt (Eq. (3)).

2. The flexibility of r_hs, taking all possible values, in contrast to Eqs. (11) and (12), where r_hs is restricted to being equal to the critical curve.

3. The simplicity and computational efficiency achieved with Eq. (13), since the elliptic integrals E,  K are eliminated.

Regarding the individual terms of Eq. (13), the second term, a_*(40 − r_hs), approximates the angular lag of a stationary emitter. As a single numerical demonstration, we find in a ray-traced hot spot with a_* = 0.5 and r_hs = 2.88 (the critical radius for this spin) that ΔPA = 198.5, near identical to δPA_0 → 1 = 198. The next term includes orbital motion through a linear dependence on Ω_coef, similarly to δPA_orb (Eq. (7)). The coefficient Ω_coef captures variation in both r_hs and a_*. We used fitting procedures to find the correct coefficients for Eq. (13).

Figure 3 illustrates the performance of this empirical relation. It compares mean ΔPA values (to avoid i dependence) obtained from simulations for a certain set of parameters with those predicted by Eq. (13). The plots have been created to show models with fixed Ω_coef (Ω_coef = 0.75 in this case) while varying (r_hs and a*). The bottom panel shows the flat difference between the two ΔPA values, which remains below 5° in all but two cases (r_hs = 4,  a_* = −0.5,   − 0.9), where it goes up to 10°.

Fig. 3. Top: Mean ΔPA values (from a full period) for different models in our library, overlaid with our empirical fitting relation. Bottom: Residuals (percentage difference between ΔPA and δPA0 → 1), where for every radius the spin values have been expanded horizontally for clarity.

In the following section we demonstrate its application through explicit examples.

4.4. Mock ΔPA observation

Equation (13) can be solved analytically for a_*, but the resulting expression is cumbersome and impractical to present here. Instead, a more convenient approach is to substitute observed values and solve for a_*. In our case, we employed a Bayesian solver using dynesty, a nested sampling tool introduced by Speagle (2020) and further developed by Koposov et al. (2023).

This choice for a Bayesian approach is motivated by two key advantages: (1) it provides a natural framework for incorporating observational uncertainties through Gaussian priors with specified standard deviations, and (2) it yields a posterior distribution for the spin parameter, inherently accounting for measurement uncertainties. By assuming a uniform prior on spin, a_* ∈ ( − 1, 1), the resulting posterior directly reflects the inferred constraints on a_*. For the likelihood calculation, we used

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{p})=-\frac{{[\mathrm{\Delta }PA-\widehat{\mathrm{\Delta }PA}(\mathit{p})]}^2}{2{\sigma }_{\mathrm{\Delta }PA}^2},\end{array}$$(15)

where p is the vector of the model parameters (P,  r_hs,  a_*), ΔPA is the observed value, while $\widehat{\mathrm{\Delta }PA}$ is the predicted value from our empirical relation in Eq. (13). The uncertainty σ_ΔPA is taken from observations.

Note that P and r_hs are observed with their own uncertainties (σ_P,  σ_rhs) but still resampled using their Gaussian distributions as prior (including the full parameter space). This is a similar approach to including all three observables in the ℒ function and only sampling a_*, but better, because it extracts more information from the system, since the posteriors can have shifted peaks from the mean of the observation. It also provides a more intuitive framework, since it samples 3 values (r_hs,  P and a_*) and extracts ΔPA that is then compared to the observation.

Given the simplicity of the task, we employed dynamic nested sampling with 1000 live points, one additional batch, and default walker and bound settings. The computation takes approximately one minute on a standard Jupyter notebook running on a typical laptop. The script for this analysis is publicly available at GitHub².

We simulated observations in which two bright spots can be identified on the screen, corresponding to the primary and secondary images of a hot spot, while following the same assumptions as our analysis modeling. A test on a more general case of non-equatorial, spiraling hot spots, is beyond the scope of the current work. In terms of accuracy, the ring diameter was measured as $\widehat{d}=51.9\pm 2\mathrm{\mu }$as in EHTC (2022a). Based on this, we assumed a comparable level of accuracy for r_hs (±0.5, M) and ΔPA (±5°). Using these uncertainties, we constructed three models, which are summarized in Table 3. The errors of ΔPA are used for the likelihood calculation of Eq. (15), while the errors of the other two parameters are used to define the prior space as a Gaussian distributions.

Before presenting the results, there is a subtle point from our empirical relation (Eq. (13)) that is worth clarifying. The way we modeled orbital velocity facilitates comparisons with numerical and analytical disk solutions; however, Ω_coef itself is not directly observable. Instead, observations provide measurements of P and r_hs, which can be related to Ω_coef and a_* through Eq. (8). This implies that we effectively have two equations with two unknowns.

To illustrate this, we present Fig. 4, where we use test 1 as a template of observation and perform two different estimations. In the left case, we estimate Ω_coef and a_* using two different fits: the “ΔPA fit”, which relies on our empirical relation (Eq. (13)), and the “P fit”, which is based on Eq. (8). The “P fit” constrains Ω_coef but leaves a_* completely unconstrained, whereas the “ΔPA fit” has no direct constraint on Ω_coef and only weakly excludes the most extreme spin values.

The bottom-middle panel on the left of Fig. 4 shows the two-dimensional posterior distribution of Ω_coef and a_*. Although each individual fit allows a_* to span nearly the full prior range, their overlap (highlighted by a black circle) demonstrates how combining the two distributions leads to a well-defined solution. This provides a clear validation of our method’s effectiveness.

Fig. 4. Bayesian fit using the mock observations, the uncertainties from Table 3, and the empirical relation of Eq. (13). Left: Fits using Eqs. (13) and (8) (period) separately. Right: Fit when combining the two.

Instead of using a two-step approach that involves overlapping fits of the two equations, we could incorporate the exploration of Ω_coef within a single algorithm, as shown on the right side of Fig. 4. In this case, Ω_coef is directly obtained from Eq. (14) and substituted into the empirical relation. Consequently, it is not treated as a free parameter in the prior and does not appear in the posteriors. Naturally, given the values of r_hs, P, and a_*, Ω_coef is fully determined. As expected from the previous exercise the spin parameter is nicely constrained at a_* = −0.16 ± 0.32 with 2σ confidence.

The results from tests 2 and 3 are presented in Fig. 5. The overall behavior remains consistent with test 1, with one notable difference: the constraining power improves when the spin values are closer to the boundaries. This is evident in the case of $a_{\ast }=0.{94}_{-0.18}^{+0.06}$, where the uncertainties are less than half those of tests 1 and 3. This behavior is expected, as the true value lies near the boundaries, but it is worth highlighting.

Fig. 5. Bayesian fit using the mock observations, the uncertainties from Table 3, and the empirical relation of Eq. (13). Note that only a* is the true constraint here; rhs and P are resampled for error propagation.

5. Discussion and conclusions

The measurement of a black hole’s spin has occupied astrophysicists for decades, even before the first black hole image from the EHT Collaboration (EHTC 2019b). With upcoming upgrades to the EHT and next-generation EHT facilities, advancements in imaging algorithms, and the development of new projects such as BHEX, novel opportunities for spin estimation have emerged (and will continue to do so). This work explores one such avenue, leveraging the interpretation of flaring events around supermassive black holes, particularly Sgr A*, as emission from a well-localized region, a hot spot. This emission is expected to be sufficiently dominant, thus making the secondary image of the hot spot detectable and opening new possibilities for probing the near-horizon space-time. Among the potential observables, we identify the difference in position angle between the primary (n = 0) and secondary (n = 1) hot spot images (ΔPA) as a key diagnostic tool.

In this work we have presented three key results that facilitate the use of this observable, under the assumptions of circular and equatorial orbits:

1. We derived the correction to the complete elliptical integrals in the analytic relation for δPA, enabling the calculation of δPA_0 → 1 with sub-degree accuracy, instead of the commonly used n → ∞ limit. Additionally, we propose a simplified approximation of the analytic solution that deviates by at most 5° at a_* → 1 for the same δPA_0 → 1 quantity.

2. We highlight a distinct effect of orbiting hot spots with regards to the system’s inclination angle (i). Specifically, when tracking the values of ΔPA over a full orbit of a moving hot spot, we find that the mean value of ΔPA remains unchanged (sub-degree for all cases apart from r_hs = 4,  a_* = −0.9, where the deviation is 3° ) for inclinations up to i ∼ 30°, with deviation of ∼5° emerging at i ∼ 60°. This finding significantly simplifies our method, as otherwise, even at i ∼ 25° – the expected inclination for Sgr A* – the fluctuations of ΔPA would require a full curve analysis and thus a more complicated comparison with models.

3. After studying a large library of more than 900 models generated using the general relativistic radiative transfer code ipole, we achieved our primary goal: constructing a simple empirical relation for estimating spin based on ΔPA observations. This method matches the analytic solution to within 5° in most cases, with a maximum deviation of 10° for r_hs ≤ 4 M. To further demonstrate the viability of this method, we conducted three tests, using mock observations with realistic uncertainties, and successfully extracted the spin in all cases, with a maximum uncertainty of ±0.3 for 2σ confidence. This is not a complete solution to the hot spot phenomenology, which can increase severely in modeling complexity, but rather a proof of concept for the usage of secondary images as space-time probes.

In the case of an observation that does not capture a full period but, for instance, 50% of it, we can still obtain valuable constraints with a more dedicated analysis. One approach is to assume an inclination angle (i ∼ 25°) and compare the ΔPA curve for several candidate values of the position angle (PA) of the black hole spin axis on the sky. This is necessary in order to model the correct phase of the curve. Independent estimates of PA include PA ∼ 130° (Ball et al. 2021; Yfantis et al. 2024a) and PA ∼ 180° (Wielgus et al. 2022b; GRAVITY Collaboration 2023; Yfantis et al. 2024b). From that we were able to extract the mean ΔPA and continue to use our method. A more robust method would be to create an adaptive Bayesian algorithm utilizing bipole, similarly to Yfantis et al. (2024a,b), and explore the full parameter space, jointly fitting for i,  PA,  r_hs,  Ω_coef, and a_*. While this approach requires significantly more computational time (up to ten days, compared to one minute for the current method), it yields a comprehensive posterior over all relevant parameters.

Looking forward, further refinement of this method could focus on expanding the parameter space to include infalling or outgoing orbits; this expansion would more readily reproduce the equatorial motions seen in general relativistic magnetohydrodynamics. Another interesting addition would be the inclusion of non-equatorial orbits; these could address a big question pertaining to the tracking of flaring structures, regarding the mechanism for hot spot creation, with off-equatorial spots often related to plasmoids rather than flux tubes. Another interesting avenue would be the creation of a machine learning algorithm that uses an expanded library that treats hot spot observational data more generally, including potentially partial sampling the time-dependent astrometry for both images (n = 0, 1) and extracting all the information for the system (i, r_hs, Ω_coef, radial velocities, etc.).

Acknowledgments

We thank the anonymous EHT internal reviewer for their comments. This publication is a part of the project Dutch Black Hole Consortium (with project number NWA 1292.19.202) of the research program of the National Science Agenda which is financed by the Dutch Research Council (NWO).

References

Appendix A: Numerical inspections

In this section we provide proof from two important tests of the simulation procedure. The first test regards the orbital velocities of the first and second images, and in particular whether the observed speeds correspond to the input on the simulation domain. The results are visible in Fig. A.2, where we show the phase of the two images, extracted from the centroid of every spot in the image domain, using the regular setup of 256x256 pixels. It is clear that the velocities in the image domain exactly follow the theoretical value, which holds even in the case of a high inclination i = 60°, where the speed changes in different parts of the orbit due to the viewing angle, but overall from start to finish the velocity is the same.

The second test, visible in Fig. A.1 is the convergence of higher order images (referring to n = 2 here) to the theoretical expectation of δPA. Additionally, it serves excellently as a demonstration of the dimness of n = 2, where the fluxes for the three spots are I₀ = 0.705,  I₁ = 0.050,  I₂ = 0.003, meaning it is virtually impossible to expect n = 2 to interfere with any observations in the near to mid future.

Fig. A.1. Orbital phase for two models with the same parameters apart from inclination: i = 0 (top) and i = 60° (bottom). The phase has been calculated using the centroids in the image domain. The black line with "x" points is the theoretical velocity used in the simulation.

Fig. A.2. Top: Snapshot from a non-moving hot spot, using a higher resolution of 1024x1024 pixels, to resolve the third image. Middle: Zoomed-in view of the third image. Bottom: ΔPA n(0 − 1) measurement from the first to the second image, and ΔPA n(1 − 2) from the second to the third. The dash-point line denotes the theoretical expectation δPA(n → ∞).

All Tables

All Figures

Fig. 1. Left: Snapshot from a hot spot movie, with resolution 256 × 256 (as in the entire library), where ΔPA is calculated from the primary to the secondary counterclockwise. Middle: Curve of ΔPA measurement with respect to time for a full period. Right: Time lag (difference in arrival times) between the first and second image. The gray vertical line in the right panels denotes the snapshot’s time in the left panels.

In the text

	Fig. 2. Distribution (violins) of ΔPA values during a full period for different inclination values. The crosses denote the average values per period and the dashed lines the value at i = 0° for a fixed a*.
In the text

	Fig. 3. Top: Mean ΔPA values (from a full period) for different models in our library, overlaid with our empirical fitting relation. Bottom: Residuals (percentage difference between ΔPA and δPA0 → 1), where for every radius the spin values have been expanded horizontally for clarity.
In the text

	Fig. 4. Bayesian fit using the mock observations, the uncertainties from Table 3, and the empirical relation of Eq. (13). Left: Fits using Eqs. (13) and (8) (period) separately. Right: Fit when combining the two.
In the text

	Fig. 5. Bayesian fit using the mock observations, the uncertainties from Table 3, and the empirical relation of Eq. (13). Note that only a* is the true constraint here; rhs and P are resampled for error propagation.
In the text

	Fig. A.1. Orbital phase for two models with the same parameters apart from inclination: i = 0 (top) and i = 60° (bottom). The phase has been calculated using the centroids in the image domain. The black line with "x" points is the theoretical velocity used in the simulation.
In the text

	Fig. A.2. Top: Snapshot from a non-moving hot spot, using a higher resolution of 1024x1024 pixels, to resolve the third image. Middle: Zoomed-in view of the third image. Bottom: ΔPA n(0 − 1) measurement from the first to the second image, and ΔPA n(1 − 2) from the second to the third. The dash-point line denotes the theoretical expectation δPA(n → ∞).
In the text