The observed radiation from hot gas accreting onto a black hole depends on both the details of the flow and the spacetime geometry. The lensing behavior of a black hole produces a distinctive pattern of autocorrelations within its photon ring that encodes its mass spin and inclination. In particular the time autocorrelation of the light curve is expected to display a series of peaks produced by light echoes of the source with each peak delayed by the characteristic time lapse inlineformulammlmath displayinlinemmlmrowmmlmimmlmimmlmrowmmlmathinlineformula between light echoes. However such peaks are absent from the light curves of observed black holes. Here we develop an analytical model for such light curves that demonstrates how even though light echoes always exist in the signal they do not produce autocorrelation peaks if the characteristic correlation timescale inlineformulammlmath displayinlinemmlmsubmmlmimmlmimmlmn0mmlmnmmlmsubmmlmathinlineformula of the source is greater than inlineformulammlmath displayinlinemmlmimmlmimmlmathinlineformula. We validate our model against simulated light curves of a stochastic accretion model ray traced with a generalrelativistic code and then fit the model to an observed light curve for Sgr inlineformulammlmath displayinlinemmlmrowmmlmsupmmlmrowmmlmi mathvariantnormalAmmlmimmlmrowmmlmrowmmlmommlmommlmrowmmlmsupmmlmrowmmlmathinlineformula. We infer that inlineformulammlmath displayinlinemmlmrowmmlmsubmmlmrowmmlmimmlmimmlmrowmmlmrowmmlmn0mmlmnmmlmrowmmlmsubmmlmogtmmlmommlmimmlmimmlmrowmmlmathinlineformula providing an explanation for the absence of light echoes in the time autocorrelations of Sgr inlineformulammlmath displayinlinemmlmsupmmlmrowmmlmi mathvariantnormalAmmlmimmlmrowmmlmommlmommlmsupmmlmathinlineformula light curves. Our results highlight the importance for black hole parameter inference of spatially resolving the photon ring via future spacebased interferometry.
['Cosmology', 'Astrophysics', 'and Gravitation', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Astrophysics of Galaxies', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Theory']
Explanation for the Absence of Secondary Peaks in Black Hole Light Curve Autocorrelations
{'No Header': "Alejandro Cárdenas-Avendaño , 1, 2, ∗ Charles Gammie , 3, † and Alexandru Lupsasca 4, ‡  \n- 1 Princeton Gravity Initiative, Princeton University, Princeton, New Jersey 08544, USA  \n2 Department of Physics, Princeton University, Princeton, New Jersey 08544, USA  \n- 3 Astronomy, Physics, NCSA, and ICASU, University of Illinois, Urbana, Illinois 61801, USA  \nThe observed radiation from hot gas accreting onto a black hole depends on both the details of the flow and the spacetime geometry. The lensing behavior of a black hole produces a distinctive pattern of autocorrelations within its photon ring that encodes its mass, spin, and inclination. In particular, the time autocorrelation of the light curve is expected to display a series of peaks produced by light echoes of the source, with each peak delayed by the characteristic time lapse τ between light echoes. However, such peaks are absent from the light curves of observed black holes. Here, we develop an analytical model for such light curves that demonstrates how, even though light echoes always exist in the signal, they do not produce autocorrelation peaks if the characteristic correlation timescale λ 0 of the source is greater than τ . We validate our model against simulated light curves of a stochastic accretion model ray traced with a general-relativistic code, and then fit the model to an observed light curve for Sgr A ∗ . We infer that λ 0 > τ , providing an explanation for the absence of light echoes in the time autocorrelations of Sgr A ∗ light curves. Our results highlight the importance for black hole parameter inference of spatially resolving the photon ring via future space-based interferometry.  \nIntroduction. -Light emitted from hot gases accreting onto black holes has been observed for decades across the electromagnetic spectrum [1-3]. This radiation depends on both the details of the astrophysical sources and the spacetime geometry around the black holes.  \nSuccessive images appearing within the photon ring are increasingly demagnified, rotated, and time-delayed. In the simplest case of an equatorial source viewed by a distant observer on the black hole spin axis, these images accumulate near the critical curve and may be labeled by the number of polar half-orbits that the corresponding photons execute around the black hole before reaching the observer. If the n th image of a point source appears at a time t n , at an angle φ n around the critical curve and at a perpendicular distance d n from it, then one can analytically prove [9] that the next image will appear at  \nMore precisely, a single source around a black hole can produce multiple images arising from photons that circumnavigate the event horizon a different number of times on their way to the observer. These mirror images are lensed into a distinctive 'photon ring' that represents the stamp imprinted on a black hole image by its strong gravity [4-7], and which tracks a 'critical curve' [8].  \nt n +1 ≈ t n + τ, φ n +1 ≈ φ n + δ, d n +1 ≈ e -γ d n , (1)  \nwhere the 'critical parameters' τ , δ , and γ -controlling the time delay, rotation and demagnification of strongly lensed images, respectively-are known functions of the black hole mass and spin [9, 10]. In particular, τ ≈ 16 M for most values of the spin, where M denotes the black hole mass and we work in geometric units with G = c = 1 .  \nBecause of the lensing behavior (1) of the black hole, the autocorrelation of its photon ring image intensity must display a distinctive multi-peaked structure, with the heights and locations of successive peaks respectively demagnified by e -γ and shifted in the spatio-temporal correlation plane (∆ t, ∆ ϕ ) by ( τ, δ ) [11] (Fig. 1 therein).  \nThe recent horizon-scale images taken by the Event Horizon Telescope (EHT) of the supermassive black holes M87 ∗ [12, 13] and Sgr A ∗ (the one at the center of our Galaxy) [14] are unable to resolve their photon rings. Vigorous efforts to extend the EHT array to space are now underway [15-17], and future observations using very-long-baseline interferometry (VLBI) to space could achieve the resolution needed to measure these rings and their predicted (ring-averaged) autocorrelations [11].  \nIn the meantime, one can already access the black hole light curves observed over many frequencies and for many sources, including active galactic nuclei, X-ray binaries, and gamma-ray bursts [18-20]. Such light curves may be regarded as (single-pixel) 'images' that are completely spatially averaged (over both radius and angle). Based on the preceding discussion, one would expect the temporal autocorrelation of many of these light curves to display multiple peaks, with each successive one demagnified by e -γ and delayed in time by the characteristic interval τ between light echoes [11, 21-25]. However, such time autocorrelations have never been detected.  \nIn particular, an analysis of a decade of 230 GHz light curves of Sgr A ∗ reported a characteristic autocorrelation timescale of 8 +3 -4 hours at 95% confidence, present down to at least a few Schwarzschild radii [26]. Assuming a mass of 4 . 3 × 10 6 M ⊙ , this value corresponds in geometric units to a timescale of 1361 +510 -680 M , which is significantly higher than the expected light echo time delay τ ≈ 16 M . An analysis of newer data collected during the 2017 EHT observation campaign of Sgr A ∗ reported a characteristic autocorrelation timescale of ∼ 1 hour, or ∼ 170 M [27].  \nDespite being significantly lower, this timescale is still much longer than the expected light echo time delay of τ ≈ 5 minutes. This raises the obvious question:", 'Where are the light echo autocorrelation peaks?': "In this letter, we revisit the theoretical expectations for the time autocorrelation of a black hole light curve and argue that the peaks caused by lensed images of the main emission are only present if the characteristic timescale λ 0 of temporal correlations in the source is much shorter than the light echo time delay τ . By contrast, if λ 0 ≳ τ , then these maxima ought to be absent, even if the lensed images are present and contribute flux to the light curve.  \nWe validate our model against simulated light curves of a stochastic accretion model that we ray trace using a general-relativistic code, and then we fit the model to an observed light curve for Sgr A ∗ . We infer that λ 0 > τ , providing an explanation for the absence of secondary peaks in the time autocorrelations of Sgr A ∗ light curves.  \nTo support this claim, we derive an analytical model for the light curve of a black hole that is surrounded by an equatorial source observed 'face-on' (that is, at a small inclination θ o from the spin axis). We then argue that the model continues to hold provided the parameter a ∗ sin θ o remains small, where a ∗ = J/M 2 ∈ [ -1 , 1] denotes the black hole spin and J its angular momentum.  \nTheoretical expectations. -Consider a polar observer ( θ o = 0 · ) of equatorial emission around the black hole. We decompose the full image into layers labeled by n ,  \nI ( t o , α, β ) = ∞ ∑ n =0 I n ( t o , α, β ) , (2)  \nwhere the n th layer corresponds to the image of the source produced by photons that travel n half-orbits around the black hole, described in Cartesian coordinates ( α, β ) on the image plane at observation time t o [8]. It is sometimes more convenient to use a polar angle φ and perpendicular distance d from the critical curve as image coordinates [9]. It follows from the lensing equations (1) that  \nI n ( t o , φ, d ) ≈ I n -1 ( t o -τ, φ -δ, e γ d ) , (3)  \nup to small corrections in 1 /n that are already negligible for n ≳ 2 . Each image layer has a flux ('light curve')  \nL n ( t o ) = ∫ I n ( t o , α, β ) d α d β. (4)  \nThe lensing relation (3) implies that  \nL n ( t o ) ≈ e -γ L n -1 ( t o -τ ) . (5)  \nHence, the total observed light curve is approximately  \nL ( t o ) = ∞ ∑ n =0 L n ( t o ) ≈ ∞ ∑ n =0 e -nγ L 0 ( t o -nτ ) . (6)  \nAs expected, it consists of a superposition of multiple copies of the light curve L 0 ( t o ) of the direct emission.  \nEach copy carries e -γ less flux and is time-delayed by τ relative to its predecessor. If the source is stationary, then the covariance C 0 (∆ t ) = ⟨ L 0 ( t o ) L 0 ( t o +∆ t ) ⟩ of the light curve for the direct emission is time-translation invariant, and the full light curve covariance at lag ∆ t is 1  \nC (∆ t ) = ⟨ L ( t o ) L ( t o +∆ t ) ⟩ (7)  \n= ∞ ∑ i =0 ∞ ∑ j =0 e -( i + j ) γ ⟨ L i ( t ) L j ( t +∆ t ) ⟩ (8)  \n≈ ∞ ∑ i =0 ∞ ∑ j =0 e -( i + j ) γ C 0 (∆ t +( i -j ) τ ) (9)  \n= 1 1 -e -2 γ ∞ ∑ m = -∞ e -| m | γ C 0 (∆ t + mτ ) , (10)  \nwhich is also stationary. This expression agrees with the conclusions arrived at by different means in Ref. [11] and describes a train of correlation peaks separated in time by τ and exponentially decreasing in height by e -γ .  \nFor moderate inclinations, the only meaningful change one can expect is that strongly lensed photons may skirt a range of bound orbits at different radii ˜ r in the 'photon shell' of trapped null geodesics [6]. The time delay τ (˜ r ) and Lyapunov exponent γ (˜ r ) are functions of this orbital radius, so one expects a smearing of Eq. (10) over a range of ˜ r , which remains quite narrow up to θ o ≲ 45 · . At ˜ r 0 , τ ≈ 16 M for all a ∗ as e -γ ranges from e -π ≈ 4% to 10% .  \nAs argued in Sec. 7 of Ref. [11], even though the lens equations (1) and (3) are modified for θ o > 0 , the formula (10) nevertheless remains exact to leading (linear) order in a ∗ sin θ o , with the first correction coming in only at subleading (quadratic) order. Thus, the only relevant time delay at small inclinations is the orbital half-period τ (˜ r 0 ) of the null geodesics trapped at the radius ˜ r 0 where bound photons have vanishing spin angular momentum. Likewise, the Lyapunov exponent γ (˜ r 0 ) governing their orbital instability fully controls the demagnification [9].  \nTherefore, we take the expression (10) (with τ and γ always evaluated near ˜ r 0 ) as our general analytical model for the covariance of the light curve. Since e -γ ≲ 10% , we expand the time autocorrelation of the light curve as  \nC (∆ t ) = C (∆ t ) C (0) ≈ C 0 (∆ t ) [ 1 -2 e -γ C 0 ( τ ) ] (11)  \n+ e -γ [ C 0 (∆ t + τ ) + C 0 (∆ t -τ )] ,  \nwhere C 0 (∆ t ) = C 0 (∆ t ) /C 0 (0) is the direct light curve time autocorrelation, and we have suppressed O ( e -nγ ) terms with n ≥ 2 , incurring only negligible errors ≲ 1% .  \nThe prediction (11) is our main theoretical result. It is analytically well-motivated and we numerically validate it in the next section. Our key point is that when the direct emission has a characteristic correlation timescale λ 0 ≳ τ , the expected peaks at ∆ t ≈ τ in Eq. (11) vanish.  \nMore precisely, Eq. (10) predicts correlation peaks at regular intervals ∆ t ∼ mτ . We now focus on the peak at ∆ t ∼ τ , expected to be produced by lensed images of the source whose photon half-orbit numbers differ by 1. Here, it is important to distinguish between two regimes.  \nIf τ ≫ λ 0 , then C 0 ( τ ) ≈ C 0 (∆ t + τ ) ≈ 0 , and thus Eq. (11) reduces to the simpler form  \nThe typical time autocorrelation C 0 (∆ t ) for a source with temporal correlations on a characteristic timescale λ 0 is a monotonically decreasing function, which starts at C 0 (0) = 1 (by definition), slowly drops until ∆ t ∼ λ 0 , and then decays exponentially for ∆ t ≳ λ 0 .  \nC (∆ t ) ≈ C 0 (∆ t ) + e -γ C 0 (∆ t -τ ) , (12)  \ndescribing a main peak at ∆ t = 0 (corresponding to the perfect autocorrelation of the signal with itself), followed by the expected secondary peak at ∆ t = τ .  \nBy contrast, if τ ≪ λ 0 , then C 0 ( τ ) ≈ 1 and moreover, C 0 (∆ t ± τ ) ≈ C 0 (∆ t ) for lags τ ∼ λ 0 , so Eq. (11) predicts  \nC (∆ t ) ∆ t ∼ τ ≈ C 0 (∆ t ) . (13)  \nHence, there should be no secondary peak in this regime.  \nState-of-the-art simulations using general-relativistic magnetohydronamics (GRMHD) typically find accretion flows that circulate around the black hole at very slightly sub-Keplerian velocities Ω ≈ ξ Ω K , with sub-Keplerianity ξ ≲ 1 [29, 30]. There is also experimental evidence for such behavior from observations with GRAVITY [31].  \nAs reported in the Supplemental Material (SM) [28], we observe these two behaviors in our numerical simulations. In the intermediate regime τ ∼ λ 0 , we find that as τ increases, the secondary peak changes from a 'bump' to an 'excess' and usually disappears well before τ ≈ λ 0 .  \nIf the characteristic timescale of correlations tracks the orbital period of the circularized flow, so that λ 0 ∼ 2 π/ Ω , then λ ≫ τ everywhere in the gas. In light of Eq. (13), this provides an explanation for the absence of secondary peaks in observed black hole light curve autocorrelations.  \nFor (i), we use inoisy [32] to simulate an equatorial source with stochastic fluctuations and AART [33] to ray trace its relativistic images. More exactly, we use inoisy to generate realizations of Gaussian random fields (GRF) with a Matérn covariance, which serve as a proxy for a hot gas that surrounds the black hole and fluctuates with a prescribed correlation structure. We then ray trace these realizations using AART , a code that exploits the integrability of light propagation in the Kerr spacetime to efficiently produce high-resolution black hole movies.  \nApplications. -We apply our model (11) for the time autocorrelation of light curves to (i) synthetic data and (ii) a light curve of Sgr A ∗ from the 2017 EHT campaign.  \nThis semi-analytical approach to producing black hole movies with an analytically known covariance function is arguably the 'optimal setup' for extracting echoes from light curves, since this method gives us complete control over all correlation scales and black hole parameters. See the SM for the details of the implementation [28].  \nWe consider a Kerr black hole with spin a ∗ = 94% observed from an inclination θ o = 20 · , and sources with different characteristic timescales λ 0 . We compare their resulting autocorrelations C (∆ t ) to the prediction (11). The blue line in Fig. 1 (top panel) shows the light curve (6) corresponding to a Keplerian flow with λ 0 = 2 π/ Ω K , i.e., proportional to r 3 / 2 s , with r s the equatorial radius in the source. For this position-dependent correlation time, the resulting GRF is inhomogeneous and anisotropic (see SM for an example snapshot [28]). The blue open crosses in Fig. 1 (bottom-left panel) plot the corresponding time autocorrelation, which does not present a secondary peak at ∆ t ∼ τ , consistent with Eq. (13) and its implications.  \nIf the source has constant correlation scales, then the light curve L 0 ( t o ) of the direct image I 0 ( t o , α, β ) is a Matérn field with d = 1 , ν = 3 / 2 and correlation length λ 0 [28]. Its autocorrelation is  \nWe stress that the absence of correlation peaks in these simulations cannot be attributed to limitations in the computation of the light curves (as we have full control of the simulation resolution and sampling rate) or in their analysis (we also applied high-pass filters and computed derivatives to search for concavity changes, to no avail).  \nC 0 (∆ t ) = ( 1 + | ∆ t | λ 0 ) exp ( -| ∆ t | λ 0 ) . (14)  \nWe emphasize that this is not a numerical fit to inoisy simulations, but rather the analytical formula describing the autocorrelation of the underlying stochastic model. This autocorrelation is also in agreement with our ray traced light curves, even when the observer is inclined and the black hole rotates rapidly, as shown in Fig. 1 (bottom-left panel) for two different values of λ 0 . The case with λ 0 = 1 . 5 M ≪ τ displays a clear secondary peak (or 'bump') while the case with λ 0 = 3 . 0 M < τ leads only to a milder 'excess' in the autocorrelation.  \nFor (ii), we compute the autocorrelation of an observed intensity light curve for Sgr A ∗ and call upon the intuition built from our model to interpret the results. Specifically, we use the April 7th, 2017 data from the 229 . 1 GHz (HI) ALMA (A1) band [27]. In Fig. 1 (top panel), we show the observed Sgr A ∗ light curve (orange points), as well as a simulated light curve computed from an M87 ∗ -like GRMHD simulation (green triangles) of a magnetically arrested disk (MAD) with r high = 40 around a black hole with a ∗ = 85% and θ o = 163 · [34]; see Refs. [34, 35] for more details. To obtain an autocorrelation from the observed Sgr A ∗ light curve, we must account for its gaps. As in Ref. [27], we compute a locally normalized discrete correlation function (LNDCF) [36, 37] and use Gaussian process regression (GPR) [38] to interpolate the data.  \nWhen the correlation timescale λ 0 varies across the source, we cannot provide an analytical expression for the resulting autocorrelation like the one in Eq. (14). We can, however, use Eq. (14) to derive an effective λ 0 by fitting it to the autocorrelation data. For the Keplerian flow with λ 0 = 2 π/ Ω K , we obtain a good fit with an effective λ 0 = 54 . 0 M , as shown with the solid blue line in Fig. 1.  \nFIG. 1. Top panel: Light curves from an observation of Sgr A ∗ with ALMA [27], from a GRMHD simulation [34], and from an inoisy simulation with λ 0 = 2 π/ Ω K ray traced with AART . We assumed a mass M = 4 . 3 × 10 6 M ⊙ for Sgr A ∗ to convert time to units of M . Bottom left: The autocorrelation of synthetic light curves produced with inoisy using different correlation times and ray traced with AART assuming a black hole with spin a ∗ = 94% observed from an inclination θ o = 20 · . For this geometry, the Lyapunov exponent γ (˜ r ) ranges over [2 . 34 , 2 . 79] and the time delay τ (˜ r ) /M ranges over [15 . 20 , 15 . 85] , where the lower and upper bounds correspond to the innermost radius ˜ r -and outermost radius ˜ r + , respectively, within the photon shell. We have used their mean values (written in the plots) for the analytical predictions presented with lines. The transparent lines represent all the possible predictions when using the whole set of values of γ and τ . The grey and pink open circles correspond to underlying correlation timescales λ 0 = 1 . 5 M and λ 0 = 3 . 0 M , respectively. The blue crosses correspond to the autocorrelation of the light curve shown in the top panel computed from an inoisy and AART simulation. Bottom right: The autocorrelations of the light curves presented in the top figure. For the observed light curve of Sgr A ∗ , we computed the autocorrelation using the LNDCF algorithm directly (orange dots), and after interpolating it with two different kernels for a Gaussian process regression (purple pluses, for Matérn, and blue crosses, for rational quadratic). The solid line corresponds to the best-fit to Eq. (14) presented in the bottom-left panel. For comparison, as in Ref. [27], we include as a dashed line an exponential decay with a 1 hour timescale. From these results, we infer that λ 0 > τ , providing an explanation for the absence of light echoes in the time autocorrelation of Sgr A ∗ light curves.  \n<!-- image -->  \nThis then allows us to apply the same procedure used for the synthetic light curves. See the SM for details of these implementations [28]. The resulting autocorrelations are similar regardless of the method used to obtain them.  \nAlthough the GRMHD and inoisy + AART light curves in Fig. 1 (top panel) are produced from very different models, their time autocorrelations (bottom-right panel) are remarkably similar to each other and also to the ones computed from the observed light curve of Sgr A ∗ .  \nConsistent with Ref. [27], Fig. 1 (bottom-right panel) shows no clear signs of a correlation peak at lag ∆ t ∼ τ that could be interpreted as an effect of lensing by Sgr A ∗ .  \nWe plot the complete autocorrelation in the SM [28]. Since C ( τ ) ≈ 1 , the plot strongly suggests that λ Sgr A ∗ 0 > τ : that is, the characteristic timescale of correlations in the plasma around Sgr A ∗ appears to exceed the time delay between light echoes. Hence, in accordance with Eq. (13), we should not expect to see secondary correlation peaks, explaining their absence from observations.  \nIn passing, we note that the Sgr A ∗ autocorrelation is well-approximated by the simulated autocorrelation (14) of a Keplerian flow with effective λ 0 = 54 . 0 M > τ (solid blue lines in the bottom panels of Fig. 1).  \nDiscussion -We have developed a simple analytical model for the autocorrelation of black hole light curves that offers insight into some of the challenges involved in separating the effects of the plasma from those of the spacetime geometry. When applied to an observed light curve of Sgr A ∗ , our model suggests that the temporal correlations inherent in its surrounding plasma suppress the autocorrelation peaks expected from lensing around the black hole, explaining the absence of such signatures.  \nWe thank Maciek Wielgus and George Wong for their valuable comments, and for providing the data for the  \nThese results indicate that the inference of black hole parameters from strong lensing effects will be difficult via light curve autocorrelations alone, and likely require future space-VLBI observations that spatially resolve the photon ring. Planning for such observations is underway.  \nSgr A ∗ and GRMHD light curves, respectively. We also thank Neal Dalal, Suvendu Giri, Lennox Keeble, Aviad Levis, Leo Stein and Sam Gralla for helpful discussions. CG and AL also thank the Aspen Center for Physics, which is supported by NSF grant PHY-2210452. ACA acknowledges support from the Simons Foundation. AL is supported by NSF grants AST-2307888 and PHY2340457. CG was supported in part by the IBM Einstein Fellow Fund at the Institute for Advanced Study. Some of the simulations presented in this work were performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University.", 'Appendix A: Variation of the photon ring critical parameters with black hole spin and inclination': 'As discussed in the main text, the presence of a black hole in the midst of an astrophysical source such as a hot radiating plasma is expected to produce light echoes in its observational appearance. 2 The time interval between light echoes can depend on the details of the source, as well as on the black hole mass, spin, and inclination. However, for strongly lensed echoes, the specific value of this time delay is fairly insensitive to these parameters. For instance, as shown in the left panel of Fig. 2, the time lapse τ incurred per half orbit is around 16 M for most values of spin and inclination [9, 21, 22], only exceeding 30 M when the black hole is near extremality and observed from high inclination. Assuming that τ takes a single value is a conservative approach. In practice, different sources produce a distribution of time delays in low-order images, broadening out the time autocorrelation and thereby obscuring light echo peaks.  \nThe effects of strong lensing on the time autocorrelation of the light curve in our model, Eq. (11), are controlled by the critical parameters γ and τ . In particular, the relativistic contributions are suppressed by factors of e -γ . The right panel of Fig. 2 illustrates that the expected suppression is approximately e -π ∼ 4% for most spins and inclinations.  \n<!-- image -->  \nFIG. 2. Left: The maximum value of the time delay τ (in units of M ) around the critical curve (or equivalently, the maximum time lapse incurred during a half-orbit within the photon shell), as a function of black hole spin and observer inclination. The time delay is τ ≈ 16 M across most of parameter space, and only grows larger for rapidly spinning black holes observed from high inclinations. The black dot represents the location of our fiducial model with spin a ∗ = 94% , inclination θ o = 20 · , and a maximal half-orbital period τ max = 15 . 85 M . Right: The minimum value of the Lyapunov exponent γ (rescaled by π ) around the critical curve (which controls the demagnification factor e -γ ), as a function of black hole spin and observer inclination. This exponent is γ ≈ π across most of parameter space, though it can decrease significantly at high spin. The white dot represents our fiducial model with spin a ∗ = 94% , inclination θ o = 20 · , and a minimum Lyapunov exponent γ min = 0 . 74 π .  \n<!-- image -->', 'Appendix B: Simulated black hole movies': "In this section, we briefly summarize the key ingredients needed to simulate the black hole movies analyzed in the main text. We first review how to ray trace black hole images with the AART code [33], and then describe our method for simulating fluctuating astrophysical sources using the inoisy code [32]. Full details are provided in Refs. [32, 33]. Ray tracing. -We ray trace black hole images-frames in a black hole movie -following the approach presented in Refs. [40-42] and implemented in AART [33]. This method assumes that the source consists of a thin disk of emitters orbiting along circular-equatorial Keplerian geodesics until they reach the radius of the innermost stable circular orbit (ISCO). Beyond this radius, they plunge into the black hole according to Cunningham's prescription [43].  \nto the observer (along n = 0 rays), and for which there is now empirical evidence [39]. The time delay between such echoes is not controlled by the period of photon orbits, but rather by the distance between the parts of the disk reflecting X-ray emission.  \nThe intensity at a particular Cartesian position ( α, β ) on the image plane of an observer is computed by (analytically) tracing the corresponding light ray backwards into the emission region. Each time a ray intersects the equatorial disk at some Boyer-Lindquist coordinates x s = ( r s , ϕ s , t s ) , the observed intensity I o ( α, β ) increases by an amount that is determined by the source emissivity I s ( x s ) . The total observed intensity at time t o is then given by [33]  \nI o ( t o , α, β ) = N ( α,β ) -1 ∑ n =0 ζ n g 3 ( r ( n ) s , α, β ) I s ( r ( n ) s , ϕ ( n ) s , t ( n ) s ) , (B1)  \nwhere x ( n ) s = x ( n ) s ( α, β ) denotes the equatorial position where the ray intersects the equatorial plane for the ( n +1) th time along its backward trajectory from image-plane position ( α, β ) at observer time t o , up to a total number N ( α, β ) along its maximal extension outside the event horizon. The redshift factor g is determined by the motion of the emitters, while ζ n is a 'fudge' factor, which we assume to be equal to 1 for n = 0 , and 1 . 5 for n ≥ 1 . We include this factor to account for the effects of the geometrical thickness of the disk [41]. In the notation of Eq. (B1), the direct image is labeled by n = 0 , while the photon ring is the collective sum of successive contributions with n ≥ 1 .  \nStochastic model for the astrophysical source. -Wenowbriefly describe our stochastic model for the source emissivity I s ( x s ) in Eq. (B1). Following Ref. [32], we numerically solve the stochastic partial differential equation  \n[1 -∇· Λ ( x s ) ∇ ] I ( x s ) = N [det Λ ] 1 4 W ( x s ) , (B2)  \nwhere x s = ( t s , x s , y s ) are Cartesian coordinates in the equatorial source (hence the subscript ' s '), N is a normalization constant, W is a standard Gaussian white noise process, and the matrix Λ controls the local covariance of the source. In general, the resulting Gaussian random field (GRF) I ( x s ) can be inhomogeneous and anisotropic, depending on the nature of the position-dependent anisotropies and correlation lengths introduced through Λ , which takes the form  \nΛ = 2 ∑ i =0 λ 2 i u i u T i . (B3)  \nHere, the vector u 0 sets the temporal correlation structure of the flow, with λ 0 its characteristic correlation time, while u 1 and u 2 set its spatial structure, with λ 1 and λ 2 the corresponding (spatial) characteristic correlation lengths.  \nGiven a realization of the GRF-namely, a numerical solution of Eq. (B2)-with periodic boundary conditions in time and Dirichlet boundary conditions in space, we create a realization of the 'standardized' field ˆ I ( x s ) defined by  \nˆ I ( x ) = I ( x ) -⟨I ( x ) ⟩ √ ⟨I 2 ( x ) ⟩ - ⟨I ( x ) ⟩ 2 , (B4)  \nwhere the angle brackets denoting an average over the source. Then given a choice of parameter σ that controls the amplitude of astrophysical fluctuations, we generate a source emissivity [33]  \nI s ( x s ) ≡ J SU ( r s ) e σ ˆ I ( x s ) -1 2 σ 2 , (B5)  \nwhere the function J SU ( r s ) serves as an envelope for the radial profile of the fluctuating source. We take it to be the Johnson's SU distribution, a function with the conveniently simple analytical form  \nJ SU ( r ; µ, ϑ, γ ) ≡ e -1 2 [ γ +arcsinh ( r -µ ϑ )] 2 √ ( r -µ ) 2 + ϑ 2 , (B6)  \nwhose three parameters µ , ϑ , and γ respectively control the location of the profile's peak, its width, and its asymmetry [42]. Given that the values of these three parameters primarily control the appearance of the time-averaged image, but not the correlations that we focus on here, we set them to the fiducial values ( µ, ϑ, γ ) = ( r -, 1 / 2 , -3 / 2) , where r -= M -√ M 2 -a 2 denotes the radius of the inner horizon. As we will show below, this particular set of parameters produces black hole images that are broadly consistent with the 2017 EHT campaign results [40].", 'Appendix C: Realizations of various Gaussian random fields': "In this section, we describe the computation of light curves corresponding to simulated observations of stochastic sources with different correlation structures. First, we consider an equatorial disk with correlation lengths that remain constant across the disk, and then turn to the case of position-dependent correlations. For that latter case, we assume that the characteristic timescale λ 0 of temporal correlations in the source follows a Keplerian scaling λ 0 ∝ 1 / Ω K , where Ω K ∝ r 3 / 2 s denotes the angular velocity of a circular-equatorial Keplerian orbit at radius r s in the source.  \n<!-- image -->  \nFIG. 3. Intensity profiles (in logarithmic scale) of single snapshots from inoisy simulations assuming constant (left) and position-dependent (right) correlation lengths. The values for the characteristic scales are shown in the images. The inset panels show ray traced images that take into account the time variability of the emission across the image of the disk (i.e., that are ray traced with 'slow light'). For both cases, we have used the envelope given in Eq. (B6) with ( µ, ϑ, γ ) = ( r -, 1 / 2 , -3 / 2) , where r -= M -√ M 2 -a 2 is the inner horizon radius, applied as presented in Eq. (B5), with a fluctuation amplitude σ = 0 . 4 .  \n<!-- image -->", '1. Constant-scale correlations': "If the correlation coefficients λ i in the matrix Λ are constant and all equal to the same correlation scale λ , then Eq. (B2) greatly simplifies, and the field it produces is a zero-mean Matérn field with covariance given by [44]  \nC ν ( x , y ) = 1 2 ν -1 Γ( ν ) ( | x -y | λ ) ν K ν ( | x -y | λ ) , (C1)  \nwhere x and y are three-dimensional position vectors in the equatorial source, Γ( x ) denotes the Gamma function, and K ν ( x ) is the modified Bessel function of the second kind and order ν . The order depends on the dimension d of the field and is given by ν = 2 -d/ 2 .  \nAs pointed out in Refs. [32, 45], one can numerically construct inhomogeneous and anisotropic GRFs by allowing Λ to vary across the domain. In particular, following Ref. [32], we take the three-dimensional unit vectors in Eq. (B3) to have the following ( t s , x s , y s ) components:  \nu 0 ( x s ) = (1 , v x ( x s ) , v y ( x s )) , (C2)  \nu 1 ( x s ) = (0 , cos θ ( x s ) , sin θ ( x s )) , (C3)  \nu 2 ( x s ) = (0 , -sin θ ( x s ) , cos θ ( x s )) . (C4)  \nThe resulting flow has a velocity characterized by the vector v = ˆ z × x . We set θ ( x s ) = arctan( y s , -x s ) + θ ∠ , where θ ∠ = 20 · sets the opening angle of spiral arms in the flow relative to the equatorial circles of constant radius r s .  \nWhile the resulting field is still a well-defined GRF arising as a solution to the stochastic partial differential equation given in Eq. (B2), it no longer has a simple analytically known covariance. Nevertheless, the resulting field can still be well-approximated by an anisotropic Matérn-like covariance function of the form  \nC ν (∆ x ) = 1 2 ν -1 Γ( ν ) s (∆ x ) ν K ν ( s (∆ x )) , (C5)  \nwith ∆ x = | x -y | and s (∆ x ) 2 = ∆ x · Λ -1 ∆ x .  \nFIG. 4. The light curve time autocorrelation resulting from different values of the temporal ( λ 0 ) and spatial ( λ 1 and λ 2 ) correlations of the GRF produced via Eq. (B5). The lines show the derived analytical prediction, Eq. (14), while the points show the numerical results. The blue crosses show the autocorrelation for a Keplerian flow and the line corresponds to a fitted value of λ 0 , which we treat as an effective correlation timescale. The spatial correlations do not affect the time autocorrelation of these light curves, as shown by the circles and crosses, which have different spatial correlations, but equal time correlations.  \n<!-- image -->  \nWith the above choices for the correlation structure of the disk, one can generate an inhomogeneous, anisotropic, time-dependent GRF I ( x s ) , and then produce an accretion-like source, provided that one applies to the GRF an envelope to control the behavior of the observed intensity. As mentioned in the previous section, we take the source to be given by Eq. (B5), choosing an amplitude σ = 0 . 4 for the fluctuations and using as an envelope the function given in Eq. (B6) with ( µ, ϑ, γ ) = ( r -, 1 / 2 , -3 / 2) , where r -= M -√ M 2 -a 2 . An example of the resulting field for λ 0 = 1 . 5 M , λ 1 = 1 . 0 M and λ 2 = 0 . 1 M is shown in the left panel of Fig. 3.  \nDespite the intricacies involved in analytically describing the behavior of inhomogeneous and anisotropic fields, if λ 0 is a constant, then one can still create a simple analytical model for the autocorrelation of the light curve L ( t ) .  \nThese simulations, and all the ones considered in this work, were run on a regular Cartesian grid ( t s , x s , y s ) of size 4096 × 1024 × 1024 . Specifically, for each of the spatial coordinates ( x s , y s ) in the equatorial disk, we uniformly placed 1024 pixels within the range [ -30 , 30] M , resulting in a resolution of about 0 . 06 M , whereas for the time coordinate t s , we placed 4096 grid points uniformly distributed within the range [0 , 5000] M , resulting in a cadence of 1 . 22 M . In Sec. C 3 below, we verify that these resolutions are sufficient to correctly compute the autocorrelations, and to ensure that we are resolving all the relevant characteristic scales.  \nAs shown in Sec. IIID of Ref. [33], one can compute several properties of these GRFs by going to momentum space via the Fourier transform  \n˜ I ( k ) = ∫ I ( x ) e -i k · x d 3 x . (C6)  \nIn the Fourier domain, one can write the autocorrelation, given by Eq. (7), as  \n〈 ˜ L ∗ ( k 0 ) ˜ L ( k ' 0 ) 〉 = 〈 ˜ I ∗ ( k ) ˜ I ( k ' ) 〉 , (C7)  \nwhere the asterisks denote complex conjugation, and we have taken the momentum vector to be k = ( k 0 , 0 , 0) . If the field I ( x ) is homogeneous, and we let x ' = x +∆ x , then its covariance in Fourier space is given by [33]  \n〈 ˜ I ∗ ( k ) ˜ I ( k ' ) 〉 = 2 πδ ( k 0 -k ' 0 ) ˜ C ν ( k 0 ) , (C8)  \nwhere ˜ C ν ( k ) is the momentum-space covariance of the field in Eq. (C1). Combining the previous two equations yields  \n〈 ˜ L ∗ ( k 0 ) ˜ L ( k ' 0 ) 〉 = 2 πδ ( k 0 -k ' 0 ) ˜ C ν ( k 0 ) . (C9)  \nTransforming this expression back to the time domain, one obtains Eq. (14) in the main text, after evaluating Eq. (C1) for only one dimension (the temporal one), with d = 1 and therefore ν = 3 / 2 , and setting λ = λ 0 .  \nMore generally, given a d -dimensional Matérn field, it is interesting to note that integrating out k < d of its dimensions results (via essentially the same argument) in a ( d -k ) -dimensional Matérn field.  \nIn Fig. 4, we compare our analytical prediction from Eq. (14) for the light curve time autocorrelation (shown with lines) against the numerically computed results obtained from GRF realizations solving Eq. (B2) (shown with points and crosses). The grey and pink symbols correspond to simulations with λ 0 = 1 . 5 M and λ 0 = 3 . 0 M , respectively. As the figure shows, despite the (unavoidable) numerical and discretization errors (and of course the fact that one is computing realizations that depend on the underlying 'activation' white noise W ), the prediction from Eq. (14) works remarkably well. In Sec. C 3 below, we check that an ensemble of realizations follows the expected statistics and that the numerical errors are under control. Lastly, the grey circles and crosses in Fig. 4 correspond to two simulations with spatial correlation length scales, λ 1 and λ 2 differing by a factor of two. Despite this change, the resulting time autocorrelation is unaffected, as expected since these simulations share the same temporal correlation timescale λ 0 .", '2. Position-dependent correlations': 'Although the choice of envelope and constant correlation structure in the previous example is enough to produce an accretion-like emission model, such a model is to simple to reproduce the features that appear in physically motivated models or in observational data. However, this model becomes capable of qualitatively reproducing several of these features, provided that one allows the characteristic correlation scales λ i to vary across the simulation domain [32, 33].  \nIn particular, if we choose the velocity of the flow to follow a Keplerian profile by setting λ 0 ∝ 1 / Ω K , and moreover allow the correlation lengths to be proportional to the source radius (so λ 1 , λ 2 ∝ r s ), then the resulting ray traced images (shown in the inset in the right panel of Fig. 3) are remarkably similar to state-of-the-art GRMHD snapshots.  \nAs mentioned in the main text, when these parameters are not constant, we cannot derive an analytical expression for the autocorrelation such as the one given in Eq. (14). We can, however, fit the model in Eq. (14) to the numerically computed autocorrelation to obtain an effective timescale λ 0 of temporal correlations. For λ 0 = 2 π/ Ω K , λ 1 = 5 r s , and λ 1 = 0 . 5 r s (a snapshot of the resulting field is shown in the right panel of Fig. 3), we find an effective value for the correlation timescale of λ 0 = 88 M , as shown in Fig. 4 with the blue lines and points. This effective timescale of the source correlations is expected to be different from the one inferred from observations of the source lensed by the black hole ( λ fitted 0 = 54 M ), whose presence relativistically broadens the autocorrelation. Although this broadening does not provide a clear signature of the light echo (nor a correlation at the relevant timescale), if one knew the time correlation of the underlying emission, then one could indirectly infer the presence of the black hole due to this shift; for further discussion of this point, see Sec. E below.  \nWe stress that our model for the autocorrelation of the light curve, Eq. (11), is valid for any autocorrelation C 0 (∆ t ) of the underlying emission. In the main text, we gave a general argument for the behavior of the model when τ ≫ λ 0 [Eq. (12)] and when τ ≪ λ 0 [Eq. (13)]. The precise behavior in the intermediate regime τ ∼ λ 0 depends on the particular functional form of the covariance of the source. For instance, if C 0 = exp( -| ∆ t | /λ 0 ) , then the light curve autocorrelation does not display a secondary maximum at ∆ t ≈ τ when τ is less than λ 0 cosh -1 ( e γ / 2) . For γ ≈ π , this means that, for this particular choice of autocorrelation, there cannot be a secondary maximum when τ ≤ 3 . 14 λ 0 .  \nAs we argued in the main text, the absence of correlation peaks in such cases cannot be attributed to limitations in the computation of the light curves. For instance, the sensitivity (or relative change) of C is exceedingly low for the autocorrelation considered above when τ/λ 0 ≪ 1 . In fact, it is exponentially insensitive, since  \nd ln C d ln τ ∣ ∣ ∣ ∣ ∆ t = τ ≈ -2 ( τ λ 0 ) 2 e -γ . (C10)', '3. Convergence of the simulations': "Our underlying emission models are sourced by an 'activation' process: the white noise W on the right-hand side of Eq. (B2). As a result, each realization of the GRF is distinct, yet still adheres to the underlying model covariance. The impact of this variation on our observable, the time autocorrelation, is shown in the left panel of Fig. 5, where we have generated five realizations for the same underlying model and resolution. As shown with dashed lines, their variance is well approximated by 1 / √ N t , where N t corresponds to the number of points in the temporal dimension.  \nFor the light curve, the spatial components of the underlying simulation are subdominant, so one only needs to resolve the temporal scale correctly. This is achieved whenever the temporal resolution is lower than the underlying characteristic timescale (i.e., dt < λ 0 ), as we show in the right panel of Fig. 5. Therein, we can see that once this condition is met, the resulting autocorrelation is unchanged, modulo the inherent stochastic variations just discussed.  \n<!-- image -->  \nFIG. 5. Left: Five different realizations of the GRF obtained from numerically solving Eq. (B2) for the same underlying correlation structure and resolution. Their variance at large time lags is well approximated by 1 / √ N t , where N t corresponds to the number of points in the temporal dimension. Right: Different temporal and spatial resolutions for a given realization of the underlying emission structure. The resolutions are quoted as N t × N x × N y . Once the temporal resolution is higher than the underlying characteristic timescale (i.e., dt < λ 0 ), the resulting autocorrelation is computed accurately.  \n<!-- image -->", 'Appendix D: Autocorrelations from observations of Sagittarius A*': "For the simulations discussed in the main text, we have tuned the parameters of the underlying inoisy simulation to mimic the high-cadence GRMHD-simulated light curve shown in the top panel of Fig. 1 in the main text. The resulting simulated light curves do not qualitatively match the observed Sgr A ∗ light curve: this is the so-called 'variability crisis' [14]. Nevertheless, thanks to the flexibility of the inoisy model, we could also tune it to match the Sgr A ∗ light curve instead by decreasing the scale σ of fluctuations. We stress, however, that σ ∼ C (0) scales out of the autocorrelation C ( δT ) = C ( δT ) /C (0) , as can be seen from the right panel of Fig. 1 in the main text.  \nFor the unevenly sampled light curve of Sgr A ∗ , we computed its time autocorrelation using the locally normalized discrete correlation function (LNDCF), defined as [36, 37]  \nLNDCF (∆ t ) = 1 N ∆ t ∑ i,j ( a i -¯ a ∆ t )] ( b j -¯ b ∆ t ) √ ( σ 2 a ∆ t -e 2 a )( σ 2 b ∆ t -e 2 b ) , (D1)  \nwhere, following closely the notation of Ref. [27], a i and b i denote the flux density measurements of the two compared light curves, and e a and e b their estimated errors. Meanwhile, ¯ a ∆ t , ¯ b ∆ t , σ a ∆ t and σ b ∆ t are respectively the flux density means and standard deviations, calculated for each lag ∆ t using only the flux density measurements that contribute to the calculation at that particular lag. Lastly, N ∆ t is the number of data pairs contributing to the lag bin ∆ t .  \nAs a check of our implementation, in Fig. 6, we reproduce the autocorrelation computed from the ALMA data for April 7th and presented in Ref. [27]. As in Ref. [27], we have also added a line representing an exponential decay with a 1 -hour timescale, and the shaded region corresponds to autocorrelation timescales between 0 . 5 and 2 hours.  \nAs discussed in the main text, we have also implemented an alternative way to calculate autocorrelations from the data: using a Gaussian process regression (GPR) [38] method to forecast and interpolate the light curve, which allows for an even sampling. In this way, we can then apply the same procedure we used in the previous sections to compute the simulated autocorrelations. This approach was pursued in Ref. [27] to model the measured light curves.  \nA GPR uses a covariance function (or kernel) to shape the prior and posterior distributions of the Gaussian process. For the observed light curve of Sgr A ∗ , we implemented two kernels (the Matérn and rational quadratic kernels) and fitted their (hyper)parameters. As seen in Eq. C1, the Matérn covariance (kernel) is characterized by ν and a lengthscale parameter λ . By contrast, the rational quadratic kernel serves as a scale mixture of several radial basis function kernels, parameterized by a length-scale λ and a scale mixture parameter α . To both models, we have added a constant kernel to scale their magnitude. The values of the parameters in these kernels are found by maximizing the log-marginal-likelihood; that is, we optimize the parameter when fitting the GPR to the data. In the top panel of Fig. 1 in the main text, we show the measured light curve, and the resulting two predictions, with their respective error for the gaps. Both of these kernels yield a similar log-likehood ( ∼ 10 4 ).  \nFIG. 6. The autocorrelation of the Sgr A ∗ light curve measured on April 7th, 2017. As in Ref. [27], we plot an exponential decay with a 1 -hour timescale, and the shaded region corresponds to autocorrelation timescales between 0 . 5 and 2 hours. This figure can be directly compared with the third panel of Fig. 9 in Ref. [27]. The structure in the resulting autocorrelation can be attributed to the limited and irregular sampling alone [27].  \n<!-- image -->", 'Appendix E: Power spectrum of light curve autocorrelations': "Here, we briefly discuss the characteristic signatures that light echoes introduce in the power spectrum (Fourier transform) of the time autocorrelation C (∆ t ) . Suppose as in Eq. (11) that the observed light curve has autocorrelation  \nC (∆ t ) = C 0 (∆ t ) [ 1 -2 e -γ C 0 ( τ ) ] + e -γ [ C 0 (∆ t + τ ) + C 0 (∆ t -τ )] + O ( e -2 γ ) , (E1)  \nwhere C 0 (∆ t ) denotes the time autocorrelation of the direct ( n = 0 ) emission, with Fourier transform  \n˜ C 0 ( k ) = ∫ C 0 (∆ t ) e -ik ∆ t d∆ t, C 0 (∆ t ) = 1 2 π ∫ ˜ C 0 ( k ) e + ik ∆ t d k. (E2)  \nThen, the Fourier transform of the full autocorrelation C (∆ t ) is  \n˜ C ( k ) = ˜ C 0 ( k ) [ 1 -2 e -γ C 0 ( τ ) ] + e -γ [ e + ikτ ˜ C 0 ( k ) + e -ikτ ˜ C 0 ( k ) ] + O ( e -2 γ ) . (E3)  \nThus, we see that the autocorrelation power spectrum in the presence of the black hole, ˜ C ( k ) , differs from the power spectrum in the absence of a black hole, ˜ C 0 ( k ) , by a multiplicative term  \n˜ C ( k ) ˜ C 0 ( k ) = 1 + 2 e -γ [cos( τk ) -C 0 ( τ )] + O ( e -2 γ ) . (E4)  \nThis calculation only accounts for the first ( n = 1 ) light echoes. We expect the n th -order light echoes to introduce additional modes oscillating with frequency nτ but exponentially suppressed by e -nγ .  \nThese characteristic oscillations in the power spectrum may provide a way to detect the presence of light echoes in black hole light curves even when they do not produce secondary peaks in the time autocorrelation. In practice, however, we expect a detection of the signature (E4) to be very challenging for a number of reasons, including:  \n- · First, as mentioned in Sec. A, different sources can produce a distribution of time delays in low-order (smalln ) images, especially at nonzero inclination, so the oscillation cos( τk ) will be smeared over a range of τ .\n- · Second, the intrinsic correlation structure of the source cannot be directly observed, so that one can directly access only C (∆ t ) or ˜ C ( k ) , but not C 0 (∆ t ) nor ˜ C 0 ( k ) . In other words, the quantity in Eq. (E4) is unobservable.\n- · Third, the source power spectrum ˜ C 0 ( k ) may itself have oscillations with frequencies near τ . Intrinsic correlations on the same timescale λ ∼ τ as the lensing echoes may very well be present in the accretion flow.  \nTo sum up, one would have to know the details of the astrophysical source, and in particular its correlation structure ˜ C 0 ( k ) , to claim that an observed oscillation with frequency τ were produced by lensing rather than the source itself.  \n- [1] M. Schmidt, Nature 197 , 1040 (1963).\n- [2] S. Bowyer, E. T. Byram, T. A. Chubb, and H. Friedman, Science 147 , 394 (1965).\n- [3] R. Narayan and E. Quataert, Nature 615 , 597 (2023), arXiv:2303.13229.\n- [4] J. P. Luminet, A&A 75 , 228 (1979).\n- [5] S. E. Gralla, D. E. Holz, and R. M. Wald, Phys. Rev. D 100 , 024018 (2019), arXiv:1906.00873.\n- [6] M. D. Johnson, A. Lupsasca, et al. , Science Advances 6 , eaaz1310 (2020), arXiv:1907.04329.\n- [7] A. Lupsasca, D. R. Mayerson, B. Ripperda, and S. Staelens, in Recent Progress on Gravity Tests: Challenges and Future Perspectives , edited by C. Bambi and A. Cárdenas-Avendaño (Springer, 2024) arXiv:2402.01290.\n- [8] J. M. Bardeen, in Black Holes (Les Astres Occlus) , edited by C. Dewitt and B. S. Dewitt (Gordon and Breach Science Publishers, 1973) pp. 215-239.\n- [9] S. E. Gralla and A. Lupsasca, Phys. Rev. D 101 , 044031 (2020), arXiv:1910.12873.\n- [10] E. Teo, General Relativity and Gravitation 35 , 1909 (2003).\n- [11] S. Hadar, M. D. Johnson, A. Lupsasca, and G. N. Wong, Phys. Rev. D 103 , 104038 (2021), arXiv:2010.03683.\n- [12] K. Akiyama et al. (Event Horizon Telescope Collaboration), ApJ 875 , L1 (2019), arXiv:1906.11238.\n- [13] K. Akiyama et al. (Event Horizon Telescope Collaboration), A&A 681 , A79 (2024).\n- [14] K. Akiyama et al. (Event Horizon Telescope Collaboration), ApJ 930 , L12 (2022).\n- [15] L. I. Gurvits et al. , Acta Astronautica 196 , 314 (2022), arXiv:2204.09144.\n- [16] P. Kurczynski et al. , in Space Telescopes and Instrumentation 2022: Optical, Infrared, and Millimeter Wave , Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12180, edited by L. E. Coyle, S. Matsuura, and M. D. Perrin (2022) p. 121800M.\n- [17] V. Kudriashov et al. , Chinese Journal of Space Science 41 , 211 (2021), arXiv:2105.06882.\n- [18] M.-H. Ulrich, L. Maraschi, and C. M. Urry, ARA&A 35 , 445 (1997).\n- [19] R. A. Remillard and J. E. McClintock, ARA&A 44 , 49 (2006), arXiv:astro-ph/0606352.\n- [20] E. Berger, ARA&A 52 , 43 (2014), arXiv:1311.2603.\n- [21] K. Fukumura and D. Kazanas, ApJ 679 , 1413 (2008), arXiv:0712.1084.\n- [22] K. Moriyama, S. Mineshige, M. Honma, and K. Akiyama, ApJ 887 , 227 (2019), arXiv:1910.10713.\n- [23] P. M. Chesler, L. Blackburn, S. S. Doeleman, M. D. Johnson, J. M. Moran, R. Narayan, and M. Wielgus, Classical and Quantum Gravity 38 , 125006 (2021), arXiv:2012.11778.\n- [24] G. N. Wong, ApJ 909 , 217 (2021), arXiv:2009.06641.\n- [25] S. Hadar, S. Harikesh, and D. Chelouche, Phys. Rev. D 107 , 124057 (2023), arXiv:2305.11247.\n- [26] J. Dexter, B. Kelly, G. C. Bower, D. P. Marrone, J. Stone, and R. Plambeck, MNRAS 442 , 2797 (2014), arXiv:1308.5968.\n- [27] M. Wielgus et al. (Event Horizon Telescope Collaboration), ApJ 930 , L19 (2022), arXiv:2207.06829.\n- [28] See Supplemental Material, which includes Refs. [39-45], at LINK, for further details on the variation of the photon ring critical parameters with black hole spin and inclination; our simulation schemes of black hole movies and Gaussian random fields; the autocorrelation from observations of Sagittarius A*; and the power spectrum of light curve autocorrelations.\n- [29] O. Porth, Y. Mizuno, Z. Younsi, and C. M. Fromm, MNRAS 502 , 2023 (2021), arXiv:2006.03658.\n- [30] N. S. Conroy, M. Bauböck, V. Dhruv, D. Lee, A. E. Broderick, C.-k. Chan, B. Georgiev, A. V. Joshi, B. Prather, and C. F. Gammie, ApJ 951 , 46 (2023), arXiv:2304.03826.\n- [31] Gravity Collaboration, A&A 684 , A200 (2024), arXiv:2401.17764.\n- [32] D. Lee and C. F. Gammie, ApJ 906 , 39 (2021), arXiv:2011.07151.\n- [33] A. Cárdenas-Avendaño, A. Lupsasca, and H. Zhu, Phys. Rev. D 107 , 043030 (2023), arXiv:2211.07469.\n- [34] G. N. Wong, L. Medeiros, A. Cárdenas-Avendaño, and J. M. Stone, 'Measuring Black Hole Light Echoes with Very Long Baseline Interferometry,' (2024), to appear.\n- [35] G. N. Wong, B. S. Prather, V. Dhruv, et al. , ApJS 259 , 64 (2022), arXiv:2202.11721.\n- [36] J. Lehar, J. N. Hewitt, D. H. Roberts, and B. F. Burke, ApJ 384 , 453 (1992).\n- [37] R. A. Edelson and J. H. Krolik, ApJ 333 , 646 (1988).\n- [38] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning (MIT Press, 2006).\n- [39] D. R. Wilkins, L. C. Gallo, E. Costantini, W. N. Brandt, and R. D. Blandford, Nature 595 , 657 (2021), arXiv:2107.13555.\n- [40] S. E. Gralla, A. Lupsasca, and D. P. Marrone, Phys. Rev. D 102 , 124004 (2020), arXiv:2008.03879.\n- [41] A. Chael, M. D. Johnson, and A. Lupsasca, ApJ 918 , 6 (2021), arXiv:2106.00683.\n- [42] H. Paugnat, A. Lupsasca, F. H. Vincent, and M. Wielgus, A&A 668 , A11 (2022), arXiv:2206.02781.\n- [43] C. T. Cunningham, ApJ 202 , 788 (1975).\n- [44] F. Lindgren, H. v. Rue, and J. Lindström, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 , 423 (2011), https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-9868.2011.00777.x.\n- [45] G.-A. Fuglstad, F. Lindgren, D. Simpson, and H. Rue, Statistica Sinica 25 , 115 (2015), arXiv:1304.6949."}