# An explanation for the absence of echoes in black hole light curve
autocorrelations

Alejandro Cárdenas-Avendaño  0000-0001-9528-1826 cardenas-
<EMAIL_ADDRESS>Princeton Gravity Initiative, Princeton University,
Princeton, New Jersey 08544, USA  Department of Physics, Princeton University,
Princeton, New Jersey 08544, USA  Charles Gammie  0000-0001-7451-8935
<EMAIL_ADDRESS>Astronomy, Physics, NCSA, and ICASU, University of
Illinois, Urbana, Illinois 61801, USA  Alexandru Lupsasca  0000-0002-1559-6965
<EMAIL_ADDRESS>Department of Physics & Astronomy,
Vanderbilt University, Nashville, Tennessee 37212, USA

###### Abstract

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 $\tau$ 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 $\lambda_{0}$ of the source
is greater than $\tau$. 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
$\lambda_{0}>\tau$, 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.

Introduction.—Light emitted from hot gases accreting onto black holes has been
observed for decades across the electromagnetic spectrum Schmidt (1963);
Bowyer  _et al._ (1965); Narayan and Quataert (2023). This radiation depends
on both the details of the astrophysical sources and the spacetime geometry
around the black holes.

More 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 Luminet (1979); Gralla  _et al._
(2019); Johnson  _et al._ (2020); Lupsasca  _et al._ (2024), and which tracks
a “critical curve” Bardeen (1973).

Successive 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^{\text{th}}$ image of a point source
appears at a time $t_{n}$, at an angle $\varphi_{n}$ around the critical curve
and at a perpendicular distance $d_{n}$ from it, then one can analytically
prove Gralla and Lupsasca (2020) that the next image will appear at

$\displaystyle t_{n+1}\approx
t_{n}+\tau,\quad\varphi_{n+1}\approx\varphi_{n}+\delta,\quad d_{n+1}\approx
e^{-\gamma}d_{n},$ (1)

where the “critical parameters” $\tau$, $\delta$, and $\gamma$—controlling the
time delay, rotation and demagnification of strongly lensed images,
respectively—are known functions of the black hole mass and spin Teo (2003);
Gralla and Lupsasca (2020). In particular, $\tau\approx 16M$ for most values
of the spin, where $M$ denotes the black hole mass and we work in geometric
units with $G=c=1$.

Because 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^{-\gamma}$ and shifted in the spatio-temporal correlation
plane $(\Delta t,\Delta\phi)$ by $(\tau,\delta)$ Hadar  _et al._ (2021) (Fig.
1 therein).

The recent horizon-scale images taken by the Event Horizon Telescope (EHT) of
the supermassive black holes M87∗ Akiyama  _et al._ (2019, 2024) and Sgr A∗
(the one at the center of our Galaxy) Akiyama  _et al._ (2022) are unable to
resolve their photon rings. Vigorous efforts to extend the EHT array to space
are now underway Gurvits  _et al._ (2022); Kurczynski  _et al._ (2022);
Kudriashov  _et al._ (2021), 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 Hadar  _et
al._ (2021).

In 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 Ulrich  _et al._ (1997); Remillard and
McClintock (2006); Berger (2014). 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^{-\gamma}$ and delayed in
time by the characteristic interval $\tau$ between light echoes Fukumura and
Kazanas (2008); Moriyama  _et al._ (2019); Wong (2021); Hadar  _et al._ (2021,
2023). However, such time autocorrelations have never been detected.

In particular, an analysis of a decade of 230 GHz light curves of Sgr A∗
reported a characteristic autocorrelation timescale of $8_{-4}^{+3}$ hours at
$95\%$ confidence, present down to at least a few Schwarzschild radii Dexter
_et al._ (2014). Assuming a mass of $4.3\times 10^{6}M_{\odot}$, this value
corresponds in geometric units to a timescale of $1361_{-680}^{+510}M$, which
is significantly higher than the expected light echo time delay $\tau\approx
16M$. An analysis of newer data collected during the 2017 EHT observation
campaign of Sgr A∗ reported a characteristic autocorrelation timescale of
$\sim 1$ hour, or $\sim 170M$ Wielgus  _et al._ (2022).

Despite being significantly lower, this timescale is still much longer than
the expected light echo time delay of $\tau\approx 5$ minutes. This raises the
obvious question:

Where are the light echoes?

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 $\lambda_{0}$ of temporal correlations in the source is much shorter
than the light echo time delay $\tau$. By contrast, if
$\lambda_{0}\gtrsim\tau$, then these peaks ought to be absent, even if the
lensed images are present and contribute flux to the light curve.

To 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 $\theta_{\rm o}$ from the spin axis). We then argue
that the model continues to hold provided the parameter $a_{*}\sin{\theta_{\rm
o}}$ remains small, where $a_{*}=J/M^{2}\in[-1,1]$ denotes the black hole spin
and $J$ its angular momentum.

We 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 $\lambda_{0}>\tau$,
thereby providing an explanation for the absence of light echoes in the time
autocorrelations of Sgr A∗ light curves.

Theoretical expectations.—Consider a polar observer ($\theta_{\rm
o}=0^{\circ}$) of equatorial emission around the black hole. We decompose the
full image into layers labeled by $n$,

$\displaystyle I(t_{\rm o},\alpha,\beta)=\sum_{n=0}^{\infty}I_{n}(t_{\rm
o},\alpha,\beta),$ (2)

where the $n^{\text{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 $(\alpha,\beta)$ on the image plane at
observation time $t_{\rm o}$ Bardeen (1973). It is sometimes more convenient
to use a polar angle $\varphi$ and perpendicular distance $d$ from the
critical curve as image coordinates Gralla and Lupsasca (2020). It follows
from the lensing equations (1) that

$\displaystyle I_{n}(t_{\rm o},\varphi,d)\approx
I_{n-1}\mathopen{}\mathclose{{}\left(t_{\rm
o}-\tau,\varphi-\delta,e^{\gamma}d}\right),$ (3)

up to small corrections in $1/n$ that are already negligible for $n\gtrsim 2$.
Each image layer has a flux (“light curve”)

$\displaystyle\mathscr{L}_{n}(t_{\rm o})=\int I_{n}(t_{\rm
o},\alpha,\beta)\mathop{}\\!\mathrm{d}\alpha\mathop{}\\!\mathrm{d}\beta.$ (4)

The lensing relation (3) implies that

$\displaystyle\mathscr{L}_{n}(t_{\rm o})\approx
e^{-\gamma}\mathscr{L}_{n-1}(t_{\rm o}-\tau).$ (5)

Hence, the total observed light curve is approximately

$\displaystyle\mathscr{L}(t_{\rm o})=\sum_{n=0}^{\infty}\mathscr{L}_{n}(t_{\rm
o})\approx\sum_{n=0}^{\infty}e^{-n\gamma}\mathscr{L}_{0}(t_{\rm o}-n\tau).$
(6)

As expected, it consists of a superposition of multiple copies of the light
curve $\mathscr{L}_{0}(t_{\rm o})$ of the direct emission. Each copy carries
$e^{-\gamma}$ less flux and is time-delayed by $\tau$ relative to its
predecessor. If the source is stationary, then the covariance $C_{0}(\Delta
t)=\mathopen{}\mathclose{{}\left\langle\mathscr{L}_{0}(t_{\rm
o})\mathscr{L}_{0}(t_{\rm o}+\Delta t)}\right\rangle$ of the light curve for
the direct emission is time-translation invariant, and the full light curve
covariance at lag $\Delta t$ is111We use
$\mathopen{}\mathclose{{}\left\langle\mathscr{L}_{i}(t)\mathscr{L}_{j}(t+\Delta
t)}\right\rangle\approx\mathopen{}\mathclose{{}\left\langle\mathscr{L}_{0}(t-i\tau)\mathscr{L}_{0}(t+\Delta
t-j\tau)}\right\rangle$ from Eq. (5), and then shift $(i,j)$ to
$(m,2s)=(i-j,i+j-|m|)$.

$\displaystyle C(\Delta t)$
$\displaystyle=\mathopen{}\mathclose{{}\left\langle\mathscr{L}(t_{\rm
o})\mathscr{L}(t_{\rm o}+\Delta t)}\right\rangle$ (7)
$\displaystyle=\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}e^{-(i+j)\gamma}\mathopen{}\mathclose{{}\left\langle\mathscr{L}_{i}(t)\mathscr{L}_{j}(t+\Delta
t)}\right\rangle$ (8)
$\displaystyle\approx\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}e^{-(i+j)\gamma}C_{0}\mathopen{}\mathclose{{}\left(\Delta
t+(i-j)\tau}\right)$ (9)
$\displaystyle=\frac{1}{1-e^{-2\gamma}}\sum_{m=-\infty}^{\infty}e^{-|m|\gamma}C_{0}\mathopen{}\mathclose{{}\left(\Delta
t+m\tau}\right),$ (10)

which is also stationary. This expression agrees with the conclusions arrived
at by different means in Ref. Hadar  _et al._ (2021) and describes a train of
correlation peaks separated in time by $\tau$ and exponentially decreasing in
height by $e^{-\gamma}$.

As argued in Sec. 7 of Ref. Hadar  _et al._ (2021), even though the lens
equations (1) and (3) are modified for $\theta_{\rm o}>0$, the formula (10)
nevertheless remains exact to leading (linear) order in $a_{*}\sin{\theta_{\rm
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 $\tau(\tilde{r}_{0})$ of the null geodesics trapped at the radius
$\tilde{r}_{0}$ where bound photons have vanishing spin angular momentum.
Likewise, the Lyapunov exponent $\gamma(\tilde{r}_{0})$ governing their
orbital instability fully controls the demagnification Gralla and Lupsasca
(2020).

For moderate inclinations, the only meaningful change one can expect is that
strongly lensed photons may skirt a range of bound orbits at different radii
$\tilde{r}$ in the “photon shell” of trapped null geodesics Johnson  _et al._
(2020). The time delay $\tau(\tilde{r})$ and Lyapunov exponent
$\gamma(\tilde{r})$ are functions of this orbital radius, so one expects a
smearing of Eq. (10) over a range of $\tilde{r}$, which remains quite narrow
up to $\theta_{\rm o}\lesssim 45^{\circ}$. At $\tilde{r}_{0}$, $\tau\approx
16M$ for all $a_{*}$ as $e^{-\gamma}$ ranges from $e^{-\pi}\approx 4\%$ to
$10\%$.

Therefore, we take the expression (10) (with $\tau$ and $\gamma$ always
evaluated near $\tilde{r}_{0}$) as our general analytical model for the
covariance of the light curve. Since $e^{-\gamma}\lesssim 10\%$, we expand the
time autocorrelation of the light curve as

$\displaystyle\mathscr{C}(\Delta t)=\frac{C(\Delta t)}{C(0)}$
$\displaystyle\approx\mathscr{C}_{0}(\Delta
t)\mathopen{}\mathclose{{}\left[1-2e^{-\gamma}\mathscr{C}_{0}(\tau)}\right]$
(11)
$\displaystyle\phantom{=}+e^{-\gamma}\mathopen{}\mathclose{{}\left[\mathscr{C}_{0}(\Delta
t+\tau)+\mathscr{C}_{0}(\Delta t-\tau)}\right],$

where $\mathscr{C}_{0}(\Delta t)=C_{0}(\Delta t)/C_{0}(0)$ is the direct light
curve time autocorrelation, and we have suppressed
$\mathcal{O}\mathopen{}\mathclose{{}\left(e^{-n\gamma}}\right)$ terms with
$n\geq 2$, incurring only negligible errors $\lesssim 1\%$.

The 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
$\lambda_{0}\gtrsim\tau$, the expected peaks at $\Delta t\approx\tau$ in Eq.
(11) vanish.

More precisely, Eq. (10) predicts correlation peaks at regular intervals
$\Delta t\sim m\tau$. We now focus on the peak at $\Delta t\sim\tau$, 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.

The typical time autocorrelation $\mathscr{C}_{0}(\Delta t)$ for a source with
temporal correlations on a characteristic timescale $\lambda_{0}$ is a
monotonically decreasing function, which starts at $\mathscr{C}_{0}(0)=1$ (by
definition), slowly drops until $\Delta t\sim\lambda_{0}$, and then decays
exponentially for $\Delta t\gtrsim\lambda_{0}$.

If $\tau\gg\lambda_{0}$, then
$\mathscr{C}_{0}(\tau)\approx\mathscr{C}_{0}(\Delta t+\tau)\approx 0$, and
thus Eq. (11) reduces to the simpler form

$\displaystyle\mathscr{C}(\Delta t)\approx\mathscr{C}_{0}(\Delta
t)+e^{-\gamma}\mathscr{C}_{0}(\Delta t-\tau),$ (12)

describing a main peak at $\Delta t=0$ (corresponding to the perfect
autocorrelation of the signal with itself), followed by the expected secondary
peak at $\Delta t=\tau$.

By contrast, if $\tau\ll\lambda_{0}$, then $\mathscr{C}_{0}(\tau)\approx 1$
and moreover, $\mathscr{C}_{0}(\Delta t\pm\tau)\approx\mathscr{C}_{0}(\Delta
t)$ for lags $\tau\sim\lambda_{0}$, so Eq. (11) predicts

$\displaystyle\mathscr{C}(\Delta t)\stackrel{{\scriptstyle\Delta
t\sim\tau}}{{\approx}}\mathscr{C}_{0}(\Delta t).$ (13)

Hence, there should be no secondary peak in this regime.

As reported in the Supplemental Material (SM), we observe these two behaviors
in our numerical simulations. In the intermediate regime
$\tau\sim\lambda_{0}$, we find that as $\tau$ increases, the secondary peak
changes from a “bump” to an “excess” and usually disappears well before
$\tau\approx\lambda_{0}$.

State-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 $\Omega\approx\xi\Omega_{\rm K}$, with
sub-Keplerianity $\xi\lesssim 1$ Porth  _et al._ (2021); Conroy  _et al._
(2023). There is also experimental evidence for such behavior from
observations with GRAVITY Gravity Collaboration (2024).

If the characteristic timescale of correlations tracks the orbital period of
the circularized flow, so that $\lambda_{0}\sim 2\pi/\Omega$, then
$\lambda\gg\tau$ 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.

Applications.—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.

For (i), we use inoisy Lee and Gammie (2021) to simulate an equatorial source
with stochastic fluctuations and AART Cárdenas-Avendaño  _et al._ (2023) 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.

This 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. We refer
readers to the SM for details of the implementation.

Figure 1: Top panel: Light curves from an observation of Sgr A∗ with ALMA
Wielgus  _et al._ (2022), from a GRMHD simulation Wong  _et al._ (2024), and
from an inoisy simulation with $\lambda_{0}=2\pi/\Omega_{\rm K}$ ray traced
with AART. We assumed a mass $M=4.3\times 10^{6}M_{\odot}$ 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 $\theta_{\rm o}=20^{\circ}$. For this geometry, the Lyapunov
exponent $\gamma(\tilde{r})$ ranges over $[2.34,2.79]$ and the time delay
$\tau(\tilde{r})/M$ ranges over $[15.20,15.85]$, where the lower and upper
bounds correspond to the innermost radius $\tilde{r}_{-}$ and outermost radius
$\tilde{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 $\gamma$ and $\tau$. The grey and pink open
circles correspond to underlying correlation timescales $\lambda_{0}=1.5M$ and
$\lambda_{0}=3.0M$, 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. Wielgus  _et al._ (2022),
we include as a dashed line an exponential decay with a 1 hour timescale. From
these results, we infer that $\lambda_{0}>\tau$, providing an explanation for
the absence of light echoes in the time autocorrelation of Sgr A∗ light
curves.

We consider a Kerr black hole with spin $a_{*}=94\%$ observed from an
inclination $\theta_{\rm o}=20^{\circ}$, and sources with different
characteristic timescales $\lambda_{0}$. We compare their resulting
autocorrelations $\mathscr{C}(\Delta t)$ to the prediction (11). The blue line
in Fig. 1 (top panel) shows the light curve (6) corresponding to a Keplerian
flow with $\lambda_{0}=2\pi/\Omega_{\rm K}$, i.e., proportional to $r_{\rm
s}^{3/2}$, with $r_{\rm s}$ the equatorial radius in the source. For this
position-dependent correlation time, the resulting GRF is inhomogeneous and
anisotropic (see the SM for an example snapshot). The blue open crosses in
Fig. 1 (bottom-left panel) plot the corresponding time autocorrelation, which
does not present a secondary peak at $\Delta t\sim\tau$, consistent with Eq.
(13) and its implications.

We 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).

When the source has constant correlation scales, the light curve
$\mathscr{L}_{0}(t_{\rm o})$ of the direct image $I_{0}(t_{\rm
o},\alpha,\beta)$ is a Matérn field with $d=1$, $\nu=3/2$ and correlation
length $\lambda_{0}$ (see the SM for a derivation). Its autocorrelation is

$\displaystyle\mathscr{C}_{0}(\Delta
t)=\mathopen{}\mathclose{{}\left(1+\frac{\mathopen{}\mathclose{{}\left|\Delta
t}\right|}{\lambda_{0}}}\right)\exp\mathopen{}\mathclose{{}\left(-\frac{\mathopen{}\mathclose{{}\left|\Delta
t}\right|}{\lambda_{0}}}\right).$ (14)

We 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 $\lambda_{0}$. The case with $\lambda_{0}=1.5M\ll\tau$ displays a
clear secondary peak (or “bump”) while the case with $\lambda_{0}=3.0M<\tau$
leads only to a milder “excess” in the autocorrelation.

When the correlation timescale $\lambda_{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
$\lambda_{0}$ by fitting it to the autocorrelation data. For the Keplerian
flow with $\lambda_{0}=2\pi/\Omega_{\rm K}$, we obtain a good fit with an
effective $\lambda_{0}=54.0M$, as shown with the solid blue line in Fig. 1.

For (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 Wielgus  _et al._ (2022). 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_{\rm high}=40$ around a black hole
with $a_{*}=85\%$ and $\theta_{\rm o}=163^{\circ}$ Wong  _et al._ (2024); see
Refs. Wong  _et al._ (2022, 2024) for more details. To obtain an
autocorrelation from the observed Sgr A∗ light curve, we must account for its
gaps. As in Ref. Wielgus  _et al._ (2022), we compute a locally normalized
discrete correlation function (LNDCF) Lehar  _et al._ (1992); Edelson and
Krolik (1988) and use Gaussian process regression (GPR) Rasmussen and Williams
(2006) to interpolate the data. This then allows us to apply the same
procedure used for the synthetic light curves. The SM provides details of
these implementations. The resulting autocorrelations are similar regardless
of the method used to obtain them.

Although 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∗.

Consistent with Ref. Wielgus  _et al._ (2022), Fig. 1 (bottom-right panel)
shows no clear signs of a correlation peak at lag $\Delta t\sim\tau$ that
could be interpreted as an effect of lensing by Sgr A∗. We plot the complete
autocorrelation in the SM. Since $\mathscr{C}(\tau)\approx 1$, the plot
strongly suggests that $\lambda^{\rm{Sgr\,A^{*}}}_{0}>\tau$: that is, the
characteristic timescale of correlations in the plasma around Sgr A∗ appears
to exceed the time delay between light echoes. Hence, as explained below Eq.
(13), we should not expect to see secondary correlation peaks, explaining
their absence from observations.

In passing, we note that the Sgr A∗ autocorrelation is well-approximated by
the simulated autocorrelation (14) of a Keplerian flow with effective
$\lambda_{0}=54.0M>\tau$ (solid blue lines in the bottom panels of Fig. 1).

Discussion—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 light echoes that are expected from lensing around the black
hole, explaining the absence of such signatures.

These results indicate that the inference of black hole parameters from strong
lensing effects will be difficult using light curve measurements alone, and
likely require future space-VLBI observations that spatially resolve the
photon ring. Planning for such observations is underway.

We thank Maciek Wielgus and George Wong for their valuable comments, and for
providing the data for the Sgr A∗ and GRMHD light curves, respectively. We
also thank Suvendu Giri, Lennox Keeble, Aviad Levis, Leo Stein and Sam Gralla
for helpful discussions. A.C.-A. acknowledges support from the Simons
Foundation. AL was supported in part by NSF grant 2307888. 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.

##

## Supplemental Material

##  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.222The light echoes that we study
herein are produced by orbiting light that is lensed into the photon ring.
Such echoes are different from those produced by X-ray “reverberations” that
are reflected from different parts of the accretion disk before traveling
directly to the observer (along $n=0$ rays), and for which there is now
empirical evidence Wilkins  _et al._ (2021). 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. 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
$\tau$ incurred per half orbit is around $16M$ for most values of spin and
inclination Fukumura and Kazanas (2008); Moriyama  _et al._ (2019); Gralla and
Lupsasca (2020), only exceeding $30M$ when the black hole is near extremality
and observed from high inclination.

The effects of strong lensing on the time autocorrelation of the light curve
in our model, Eq. (11), are controlled by the critical parameters $\gamma$ and
$\tau$. In particular, the relativistic contributions are suppressed by
factors of $e^{-\gamma}$. The right panel of Fig. 2 illustrates that the
expected suppression is approximately $e^{-\pi}\sim 4\%$ for most spins and
inclinations.

Figure 2: Left: The maximum value of the time delay $\tau$ (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 $\tau\approx 16M$ 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 $\theta_{\rm
o}=20^{\circ}$, and a maximal half-orbital period $\tau_{\rm max}=15.85M$.
Right: The minimum value of the Lyapunov exponent $\gamma$ (rescaled by $\pi$)
around the critical curve (which controls the demagnification factor
$e^{-\gamma}$), as a function of black hole spin and observer inclination.
This exponent is $\gamma\approx\pi$ 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 $\theta_{\rm o}=20^{\circ}$, and a
minimum Lyapunov exponent $\gamma_{\rm min}=0.74\pi$.

##  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 Cárdenas-Avendaño  _et al._ (2023),
and then describe our method for simulating fluctuating astrophysical sources
using the inoisy code Lee and Gammie (2021). Full details are provided in
Refs. Lee and Gammie (2021); Cárdenas-Avendaño  _et al._ (2023).

Ray tracing.—We ray trace black hole images—frames in a black hole
movie—following the approach presented in Refs. Gralla  _et al._ (2020); Chael
_et al._ (2021); Paugnat  _et al._ (2022) and implemented in AART Cárdenas-
Avendaño  _et al._ (2023). 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 Cunningham (1975).

The intensity at a particular Cartesian position $(\alpha,\beta)$ 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 $\mathbf{x}_{\rm
s}=\mathopen{}\mathclose{{}\left(r_{\rm s},\phi_{\rm s},t_{\rm s}}\right)$,
the observed intensity $I_{\rm o}(\alpha,\beta)$ increases by an amount that
is determined by the source emissivity $I_{\rm s}(\mathbf{x}_{\rm s})$. The
total observed intensity at time $t_{\rm o}$ is then given by Cárdenas-
Avendaño  _et al._ (2023)

$\displaystyle I_{\rm o}(t_{\rm
o},\alpha,\beta)=\sum_{n=0}^{N(\alpha,\beta)-1}\zeta_{n}g^{3}\mathopen{}\mathclose{{}\left(r_{\rm
s}^{(n)},\alpha,\beta}\right)I_{\rm s}\mathopen{}\mathclose{{}\left(r_{\rm
s}^{(n)},\phi_{\rm s}^{(n)},t_{\rm s}^{(n)}}\right),$ (15)

where $\mathbf{x}_{\rm s}^{(n)}=\mathbf{x}_{\rm s}^{(n)}(\alpha,\beta)$
denotes the equatorial position where the ray intersects the equatorial plane
for the $(n+1)^{\text{th}}$ time along its backward trajectory from image-
plane position $(\alpha,\beta)$ at observer time $t_{\rm o}$, up to a total
number $N(\alpha,\beta)$ along its maximal extension outside the event
horizon. The redshift factor $g$ is determined by the motion of the emitters,
while $\zeta_{n}$ is a “fudge” factor, which we assume to be equal to $1$ for
$n=0$, and $1.5$ for $n\geq 1$. We include this factor to account for the
effects of the geometrical thickness of the disk Chael  _et al._ (2021). In
the notation of Eq. (15), the direct image is labeled by $n=0$, while the
photon ring is the collective sum of successive contributions with $n\geq 1$.

Stochastic model for the astrophysical source.—We now briefly describe our
stochastic model for the source emissivity $I_{\rm s}(\mathbf{x}_{\rm s})$ in
Eq. (15). Following Ref. Lee and Gammie (2021), we numerically solve the
stochastic partial differential equation

$\displaystyle\mathopen{}\mathclose{{}\left[1-\nabla\cdot\mathbf{\Lambda}\mathopen{}\mathclose{{}\left(\mathbf{x}_{\rm
s}}\right)\nabla}\right]\mathcal{I}(\mathbf{x}_{\rm
s})=\mathcal{N}\mathopen{}\mathclose{{}\left[\det{\mathbf{\Lambda}}}\right]^{\frac{1}{4}}\mathcal{W}(\mathbf{x}_{\rm
s}),$ (16)

where $\mathbf{x}_{\rm s}=(t_{\rm s},x_{\rm s},y_{\rm s})$ are Cartesian
coordinates in the equatorial source (hence the subscript ‘$s$’),
$\mathcal{N}$ is a normalization constant, $\mathcal{W}$ is a standard
Gaussian white noise process, and the matrix $\mathbf{\Lambda}$ controls the
local covariance of the source. In general, the resulting Gaussian random
field (GRF) $\mathcal{I}\mathopen{}\mathclose{{}\left(\mathbf{x}_{\rm
s}}\right)$ can be inhomogeneous and anisotropic, depending on the nature of
the position-dependent anisotropies and correlation lengths introduced through
$\mathbf{\Lambda}$, which takes the form

$\displaystyle\mathbf{\Lambda}=\sum_{i=0}^{2}\lambda_{i}^{2}\mathbf{u}_{i}\mathbf{u}_{i}^{T}.$
(17)

Here, the vector $\mathbf{u}_{0}$ sets the temporal correlation structure of
the flow, with $\lambda_{0}$ its characteristic correlation time, while
$\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ set its spatial structure, with
$\lambda_{1}$ and $\lambda_{2}$ the corresponding (spatial) characteristic
correlation lengths.

Given a realization of the GRF—namely, a numerical solution of Eq. (16)—with
periodic boundary conditions in time and Dirichlet boundary conditions in
space, we create a realization of the “standardized” field
$\hat{\mathcal{I}}(\mathbf{x}_{\rm s})$ defined by

$\displaystyle\hat{\mathcal{I}}(\mathbf{x})=\frac{\mathcal{I}(\mathbf{x})-\mathopen{}\mathclose{{}\left\langle\mathcal{I}(\mathbf{x})}\right\rangle}{\sqrt{\mathopen{}\mathclose{{}\left\langle\mathcal{I}^{2}(\mathbf{x})}\right\rangle-\mathopen{}\mathclose{{}\left\langle\mathcal{I}(\mathbf{x})}\right\rangle^{2}}},$
(18)

where the angle brackets denoting an average over the source. Then given a
choice of parameter $\sigma$ that controls the amplitude of astrophysical
fluctuations, we generate a source emissivity Cárdenas-Avendaño  _et al._
(2023)

$\displaystyle I_{\rm s}(\mathbf{x}_{\rm s})\equiv J_{\rm{SU}}(r_{\rm
s})e^{\sigma\hat{\mathcal{I}}(\mathbf{x}_{\rm s})-\frac{1}{2}\sigma^{2}},$
(19)

where the function $J_{\rm{SU}}(r_{\rm 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

$\displaystyle
J_{\rm{SU}}(r;\mu,\vartheta,\gamma)\equiv\frac{e^{-\frac{1}{2}\mathopen{}\mathclose{{}\left[\gamma+\operatorname{arcsinh}\mathopen{}\mathclose{{}\left(\frac{r-\mu}{\vartheta}}\right)}\right]^{2}}}{\sqrt{\mathopen{}\mathclose{{}\left(r-\mu}\right)^{2}+\vartheta^{2}}},$
(20)

whose three parameters $\mu$, $\vartheta$, and $\gamma$ respectively control
the location of the profile’s peak, its width, and its asymmetry Paugnat  _et
al._ (2022). 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
$(\mu,\vartheta,\gamma)=(r_{-},1/2,-3/2)$, where $r_{-}=M-\sqrt{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 Gralla  _et al._ (2020).

##  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 $\lambda_{0}$ of temporal correlations in the source
follows a Keplerian scaling $\lambda_{0}\propto 1/\Omega_{\rm K}$, where
$\Omega_{\rm K}\propto r_{\rm s}^{3/2}$ denotes the angular velocity of a
circular-equatorial Keplerian orbit at radius $r_{\rm s}$ in the source.

Figure 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. (20) with $(\mu,\vartheta,\gamma)=(r_{-},1/2,-3/2)$, where
$r_{-}=M-\sqrt{M^{2}-a^{2}}$ is the inner horizon radius, applied as presented
in Eq. (19), with a fluctuation amplitude $\sigma=0.4$.

###  C.1 Constant-scale correlations

If the correlation coefficients $\lambda_{i}$ in the matrix $\mathbf{\Lambda}$
are constant and all equal to the same correlation scale $\lambda$, then Eq.
(16) greatly simplifies, and the field it produces is a zero-mean Matérn field
with covariance given by Lindgren  _et al._ (2011)

$\displaystyle
C_{\nu}(\mathbf{x},\mathbf{y})=\frac{1}{2^{\nu-1}\Gamma(\nu)}\mathopen{}\mathclose{{}\left(\frac{|\mathbf{x}-\mathbf{y}|}{\lambda}}\right)^{\nu}K_{\nu}\mathopen{}\mathclose{{}\left(\frac{|\mathbf{x}-\mathbf{y}|}{\lambda}}\right),$
(21)

where $\mathbf{x}$ and $\mathbf{y}$ are three-dimensional position vectors in
the equatorial source, $\Gamma(x)$ denotes the Gamma function, and
$K_{\nu}(x)$ is the modified Bessel function of the second kind and order
$\nu$. The order depends on the dimension $d$ of the field and is given by
$\nu=2-d/2$.

As pointed out in Refs. Fuglstad  _et al._ (2015); Lee and Gammie (2021), one
can numerically construct inhomogeneous and anisotropic GRFs by allowing
$\mathbf{\Lambda}$ to vary across the domain. In particular, following Ref.
Lee and Gammie (2021), we take the three-dimensional unit vectors in Eq. (17)
to have the following $(t_{\rm s},x_{\rm s},y_{\rm s})$ components:

$\displaystyle\mathbf{u}_{0}(\mathbf{x}_{\rm s})$
$\displaystyle=\mathopen{}\mathclose{{}\left(1,v_{x}(\mathbf{x}_{\rm
s}),v_{y}(\mathbf{x}_{\rm s})}\right),$ (22)
$\displaystyle\mathbf{u}_{1}(\mathbf{x}_{\rm s})$
$\displaystyle=\mathopen{}\mathclose{{}\left(0,\cos{\theta(\mathbf{x}_{\rm
s})},\sin{\theta(\mathbf{x}_{\rm s})}}\right),$ (23)
$\displaystyle\mathbf{u}_{2}(\mathbf{x}_{\rm s})$
$\displaystyle=\mathopen{}\mathclose{{}\left(0,-\sin{\theta(\mathbf{x}_{\rm
s})},\cos{\theta(\mathbf{x}_{\rm s})}}\right).$ (24)

The resulting flow has a velocity characterized by the vector
$\mathbf{v}=\hat{\mathbf{z}}\times\mathbf{x}$. We set $\theta(\mathbf{x}_{\rm
s})=\arctan(y_{\rm s},-x_{\rm s})+\theta_{\angle}$, where
$\theta_{\angle}=20^{\circ}$ sets the opening angle of spiral arms in the flow
relative to the equatorial circles of constant radius $r_{\rm s}$.

While the resulting field is still a well-defined GRF arising as a solution to
the stochastic partial differential equation given in Eq. (16), 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

$\displaystyle
C_{\nu}(\Delta\mathbf{x})=\frac{1}{2^{\nu-1}\Gamma(\nu)}s\mathopen{}\mathclose{{}\left(\Delta\mathbf{x}}\right)^{\nu}K_{\nu}\mathopen{}\mathclose{{}\left(s\mathopen{}\mathclose{{}\left(\Delta\mathbf{x}}\right)}\right),$
(25)

with $\Delta\mathbf{x}=|\mathbf{x}-\mathbf{y}|$ and
$s(\Delta\mathbf{x})^{2}=\Delta\mathbf{x}\cdot\mathbf{\Lambda}^{-1}\Delta\mathbf{x}$.

With the above choices for the correlation structure of the disk, one can
generate an inhomogeneous, anisotropic, time-dependent GRF
$\mathcal{I}(\mathbf{x}_{\rm 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. (19), choosing an amplitude $\sigma=0.4$ for the
fluctuations and using as an envelope the function given in Eq. (20) with
$(\mu,\vartheta,\gamma)=(r_{-},1/2,-3/2)$, where $r_{-}=M-\sqrt{M^{2}-a^{2}}$.
An example of the resulting field for $\lambda_{0}=1.5M$, $\lambda_{1}=1.0M$
and $\lambda_{2}=0.1M$ is shown in the left panel of Fig. 3.

These simulations, and all the ones considered in this work, were run on a
regular Cartesian grid $(t_{\rm s},x_{\rm s},y_{\rm s})$ of size $4096\times
1024\times 1024$. Specifically, for each of the spatial coordinates $(x_{\rm
s},y_{\rm s})$ in the equatorial disk, we uniformly placed $1024$ pixels
within the range $[-30,30]M$, resulting in a resolution of about $0.06M$,
whereas for the time coordinate $t_{\rm s}$, we placed $4096$ grid points
uniformly distributed within the range $[0,5000]M$, resulting in a cadence of
$1.22M$. 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.

Despite the intricacies involved in analytically describing the behavior of
inhomogeneous and anisotropic fields, if $\lambda_{0}$ is a constant, then one
can still create a simple analytical model for the autocorrelation of the
light curve $\mathscr{L}(t)$.

As shown in Sec. IIID of Ref. Cárdenas-Avendaño  _et al._ (2023), one can
compute several properties of these GRFs by going to momentum space via the
Fourier transform

$\displaystyle\tilde{\mathcal{I}}(\mathbf{k})=\int\mathcal{I}(\mathbf{x})e^{-i\mathbf{k}\cdot\mathbf{x}}\mathop{}\\!\mathrm{d}^{3}\mathbf{x}.$
(26)

In the Fourier domain, one can write the autocorrelation, given by Eq. (7), as

$\displaystyle\mathopen{}\mathclose{{}\left\langle\tilde{\mathscr{L}}^{*}(k_{0})\tilde{\mathscr{L}}\mathopen{}\mathclose{{}\left(k_{0}^{\prime}}\right)}\right\rangle=\mathopen{}\mathclose{{}\left\langle\tilde{\mathcal{I}}^{*}(\mathbf{k})\tilde{\mathcal{I}}\mathopen{}\mathclose{{}\left(\mathbf{k}^{\prime}}\right)}\right\rangle,$
(27)

where the asterisks denote complex conjugation, and we have taken the momentum
vector to be $\mathbf{k}=(k_{0},0,0)$. If the field
$\mathcal{I}\mathopen{}\mathclose{{}\left(\mathbf{x}}\right)$ is homogeneous,
and we let $\mathbf{x}^{\prime}=\mathbf{x}+\Delta\mathbf{x}$, then its
covariance in Fourier space is given by Cárdenas-Avendaño  _et al._ (2023)

$\displaystyle\mathopen{}\mathclose{{}\left\langle\tilde{\mathcal{I}}^{*}(\mathbf{k})\tilde{\mathcal{I}}\mathopen{}\mathclose{{}\left(\mathbf{k}^{\prime}}\right)}\right\rangle=2\pi\delta\mathopen{}\mathclose{{}\left(k_{0}-k_{0}^{\prime}}\right)\tilde{C}_{\nu}\mathopen{}\mathclose{{}\left(k_{0}}\right),$
(28)

where $\tilde{C}_{\nu}\mathopen{}\mathclose{{}\left(\mathbf{k}}\right)$ is the
momentum-space covariance of the field in Eq. (21). Combining the previous two
equations yields

$\displaystyle\mathopen{}\mathclose{{}\left\langle\tilde{\mathscr{L}}^{*}(k_{0})\tilde{\mathscr{L}}\mathopen{}\mathclose{{}\left(k_{0}^{\prime}}\right)}\right\rangle=2\pi\delta\mathopen{}\mathclose{{}\left(k_{0}-k_{0}^{\prime}}\right)\tilde{C}_{\nu}\mathopen{}\mathclose{{}\left(k_{0}}\right).$
(29)

Transforming this expression back to the time domain, one obtains Eq. (14) in
the main text, after evaluating Eq. (21) for only one dimension (the temporal
one), with $d=1$ and therefore $\nu=3/2$, and setting $\lambda=\lambda_{0}$.

Figure 4: The light curve time autocorrelation resulting from different values
of the temporal ($\lambda_{0}$) and spatial ($\lambda_{1}$ and $\lambda_{2}$)
correlations of the GRF produced via Eq. (19). 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 $\lambda_{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.

More 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.

In 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. (16) (shown with points and
crosses). The grey and pink symbols correspond to simulations with
$\lambda_{0}=1.5M$ and $\lambda_{0}=3.0M$, 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 $\mathcal{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, $\lambda_{1}$ and
$\lambda_{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 $\lambda_{0}$.

###  C.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 $\lambda_{i}$ to vary across the
simulation domain Lee and Gammie (2021); Cárdenas-Avendaño  _et al._ (2023).

In particular, if we choose the velocity of the flow to follow a Keplerian
profile by setting $\lambda_{0}\propto 1/\Omega_{\rm K}$, and moreover allow
the correlation lengths to be proportional to the source radius (so
$\lambda_{1},\lambda_{2}\propto r_{\rm 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.

As 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
$\lambda_{0}$ of temporal correlations. For $\lambda_{0}=2\pi/\Omega_{\rm K}$,
$\lambda_{1}=5r_{\rm s}$, and $\lambda_{1}=0.5r_{\rm 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 $\lambda_{0}=88M$, 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 ($\lambda_{0}^{\rm
fitted}=54M$), 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.

We stress that our model for the autocorrelation of the light curve, Eq. (11),
is valid for any autocorrelation $\mathscr{C}_{0}(\Delta t)$ of the underlying
emission. In the main text, we gave a general argument for the behavior of the
model when $\tau\gg\lambda_{0}$ [Eq. (12)] and when $\tau\ll\lambda_{0}$ [Eq.
(13)]. The precise behavior in the intermediate regime $\tau\sim\lambda_{0}$
depends on the particular functional form of the covariance of the source. For
instance, if $\mathscr{C}_{0}=\exp\mathopen{}\mathclose{{}\left(-|\Delta
t|/\lambda_{0}}\right)$, then the light curve autocorrelation does not display
a secondary maximum at $\Delta t\approx\tau$ when $\tau$ is less than
$\lambda_{0}\cosh^{-1}(e^{\gamma}/2)$. For $\gamma\approx\pi$, this means
that, for this particular choice of autocorrelation, there cannot be a
secondary maximum when $\tau\leq 3.14\lambda_{0}$.

As 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 $\mathscr{C}$ is
exceedingly low for the autocorrelation considered above when
$\tau/\lambda_{0}\ll 1$. In fact, it is exponentially insensitive, since

$\displaystyle\mathopen{}\mathclose{{}\left.\frac{d\ln\mathscr{C}}{d\ln\tau}}\right|_{\Delta
t=\tau}\approx-2\mathopen{}\mathclose{{}\left(\frac{\tau}{\lambda_{0}}}\right)^{2}e^{-\gamma}.$
(30)

###  C.3 Convergence of the simulations

Our underlying emission models are sourced by an “activation” process: the
white noise $\mathcal{W}$ on the right-hand side of Eq. (16). 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/\sqrt{N_{t}}$, where
$N_{t}$ corresponds to the number of points in the temporal dimension.

Figure 5: Left: Five different realizations of the GRF obtained from
numerically solving Eq. (16) for the same underlying correlation structure and
resolution. Their variance at large time lags is well approximated by
$1/\sqrt{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}\times N_{x}\times N_{y}$. Once the temporal resolution is
higher than the underlying characteristic timescale (i.e., $dt<\lambda_{0}$),
the resulting autocorrelation is computed accurately.

For 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<\lambda_{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.

##  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” Akiyama  _et al._ (2022).
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
$\sigma$ of fluctuations. We stress, however, that $\sigma\sim C(0)$ scales
out of the autocorrelation $\mathscr{C}(\delta T)=C(\delta T)/C(0)$, as can be
seen from the right panel of Fig. 1 in the main text.

For the unevenly sampled light curve of Sgr A∗, we computed its time
autocorrelation using the locally normalized discrete correlation function
(LNDCF), defined as Lehar  _et al._ (1992); Edelson and Krolik (1988)

$\displaystyle\textrm{LNDCF}(\Delta t)=\frac{1}{N_{\Delta
t}}\sum_{i,j}\frac{\mathopen{}\mathclose{{}\left(a_{i}-\bar{a}_{\Delta
t}}\right)]\mathopen{}\mathclose{{}\left(b_{j}-\bar{b}_{\Delta
t}}\right)}{\sqrt{\mathopen{}\mathclose{{}\left(\sigma_{a\Delta
t}^{2}-e_{a}^{2}}\right)\mathopen{}\mathclose{{}\left(\sigma_{b\Delta
t}^{2}-e_{b}^{2}}\right)}},$ (31)

where, following closely the notation of Ref. Wielgus  _et al._ (2022),
$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,
$\bar{a}_{\Delta t}$, $\bar{b}_{\Delta t}$, $\sigma_{a\Delta t}$ and
$\sigma_{b\Delta t}$ are respectively the flux density means and standard
deviations, calculated for each lag $\Delta t$ using only the flux density
measurements that contribute to the calculation at that particular lag.
Lastly, $N_{\Delta t}$ is the number of data pairs contributing to the lag bin
$\Delta t$.

As a check of our implementation, in Fig. 6, we reproduce the autocorrelation
computed from the ALMA data for April 7th and presented in Ref. Wielgus  _et
al._ (2022). As in Ref. Wielgus  _et al._ (2022), 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.

As discussed in the main text, we have also implemented an alternative way to
calculate autocorrelations from the data: using a Gaussian process regression
(GPR) Rasmussen and Williams (2006) 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. Wielgus  _et al._ (2022)
to model the measured light curves.

A 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. 21, the Matérn covariance
(kernel) is characterized by $\nu$ and a length-scale parameter $\lambda$. By
contrast, the rational quadratic kernel serves as a scale mixture of several
radial basis function kernels, parameterized by a length-scale $\lambda$ and a
scale mixture parameter $\alpha$. 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 ($\sim 10^{4}$).

Figure 6: The autocorrelation of the Sgr A∗ light curve measured on April 7th,
2017. As in Ref. Wielgus  _et al._ (2022), 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. Wielgus  _et al._ (2022). The structure
in the resulting autocorrelation can be attributed to the limited and
irregular sampling alone Wielgus  _et al._ (2022).

## References

  * Schmidt (1963) M. Schmidt, Nature 197, 1040 (1963).
  * Bowyer  _et al._ (1965) S. Bowyer, E. T. Byram, T. A. Chubb, and H. Friedman, Science 147, 394 (1965).
  * Narayan and Quataert (2023) R. Narayan and E. Quataert, Nature 615, 597 (2023), arXiv:2303.13229 [astro-ph.HE] .
  * Luminet (1979) J. P. Luminet, A&A 75, 228 (1979).
  * Gralla  _et al._ (2019) S. E. Gralla, D. E. Holz, and R. M. Wald, Phys. Rev. D 100, 024018 (2019), arXiv:1906.00873 [astro-ph.HE] .
  * Johnson  _et al._ (2020) M. D. Johnson, A. Lupsasca, _et al._ , Science Advances 6, eaaz1310 (2020), arXiv:1907.04329 [astro-ph.IM] .
  * Lupsasca  _et al._ (2024) 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 [gr-qc] .
  * Bardeen (1973) 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.
  * Gralla and Lupsasca (2020) S. E. Gralla and A. Lupsasca, Phys. Rev. D 101, 044031 (2020), arXiv:1910.12873 [gr-qc] .
  * Teo (2003) E. Teo, General Relativity and Gravitation 35, 1909 (2003).
  * Hadar  _et al._ (2021) S. Hadar, M. D. Johnson, A. Lupsasca, and G. N. Wong, Phys. Rev. D 103, 104038 (2021), arXiv:2010.03683 [gr-qc] .
  * Akiyama  _et al._ (2019) K. Akiyama _et al._ (Event Horizon Telescope Collaboration), ApJ 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Akiyama  _et al._ (2024) K. Akiyama _et al._ (Event Horizon Telescope Collaboration), A&A 681, A79 (2024).
  * Akiyama  _et al._ (2022) K. Akiyama _et al._ (Event Horizon Telescope Collaboration), ApJ 930, L12 (2022).
  * Gurvits  _et al._ (2022) L. I. Gurvits _et al._ , Acta Astronautica 196, 314 (2022), arXiv:2204.09144 [astro-ph.IM] .
  * Kurczynski  _et al._ (2022) 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.
  * Kudriashov  _et al._ (2021) V. Kudriashov _et al._ , Chinese Journal of Space Science 41, 211 (2021), arXiv:2105.06882 [astro-ph.IM] .
  * Ulrich  _et al._ (1997) M.-H. Ulrich, L. Maraschi, and C. M. Urry, ARA&A 35, 445 (1997).
  * Remillard and McClintock (2006) R. A. Remillard and J. E. McClintock, ARA&A 44, 49 (2006), arXiv:astro-ph/0606352 [astro-ph] .
  * Berger (2014) E. Berger, ARA&A 52, 43 (2014), arXiv:1311.2603 [astro-ph.HE] .
  * Fukumura and Kazanas (2008) K. Fukumura and D. Kazanas, ApJ 679, 1413 (2008), arXiv:0712.1084 [astro-ph] .
  * Moriyama  _et al._ (2019) K. Moriyama, S. Mineshige, M. Honma, and K. Akiyama, ApJ 887, 227 (2019), arXiv:1910.10713 [astro-ph.HE] .
  * Wong (2021) G. N. Wong, ApJ 909, 217 (2021), arXiv:2009.06641 [astro-ph.HE] .
  * Hadar  _et al._ (2023) S. Hadar, S. Harikesh, and D. Chelouche, Phys. Rev. D 107, 124057 (2023), arXiv:2305.11247 [gr-qc] .
  * Dexter  _et al._ (2014) J. Dexter, B. Kelly, G. C. Bower, D. P. Marrone, J. Stone, and R. Plambeck, MNRAS 442, 2797 (2014), arXiv:1308.5968 [astro-ph.HE] .
  * Wielgus  _et al._ (2022) M. Wielgus _et al._ (Event Horizon Telescope Collaboration), ApJ 930, L19 (2022), arXiv:2207.06829 [astro-ph.HE] .
  * Porth  _et al._ (2021) O. Porth, Y. Mizuno, Z. Younsi, and C. M. Fromm, MNRAS 502, 2023 (2021), arXiv:2006.03658 [astro-ph.HE] .
  * Conroy  _et al._ (2023) 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 [astro-ph.HE] .
  * Gravity Collaboration (2024) Gravity Collaboration, A&A 684, A200 (2024), arXiv:2401.17764 [astro-ph.SR] .
  * Lee and Gammie (2021) D. Lee and C. F. Gammie, ApJ 906, 39 (2021), arXiv:2011.07151 [astro-ph.IM] .
  * Cárdenas-Avendaño  _et al._ (2023) A. Cárdenas-Avendaño, A. Lupsasca, and H. Zhu, Phys. Rev. D 107, 043030 (2023), arXiv:2211.07469 [gr-qc] .
  * Wong  _et al._ (2024) G. N. Wong, L. Medeiros, A. Cárdenas-Avendaño, and J. M. Stone, “Measuring black hole mass from light echoes with very-long-baseline interfermetry,” (2024), to appear.
  * Wong  _et al._ (2022) G. N. Wong, B. S. Prather, V. Dhruv, _et al._ , ApJS 259, 64 (2022), arXiv:2202.11721 [astro-ph.HE] .
  * Lehar  _et al._ (1992) J. Lehar, J. N. Hewitt, D. H. Roberts, and B. F. Burke, ApJ 384, 453 (1992).
  * Edelson and Krolik (1988) R. A. Edelson and J. H. Krolik, ApJ 333, 646 (1988).
  * Rasmussen and Williams (2006) C. E. Rasmussen and C. K. I. Williams,  _Gaussian Processes for Machine Learning_ (MIT Press, 2006).
  * Wilkins  _et al._ (2021) D. R. Wilkins, L. C. Gallo, E. Costantini, W. N. Brandt, and R. D. Blandford, Nature 595, 657 (2021), arXiv:2107.13555 [astro-ph.HE] .
  * Gralla  _et al._ (2020) S. E. Gralla, A. Lupsasca, and D. P. Marrone, Phys. Rev. D 102, 124004 (2020), arXiv:2008.03879 [gr-qc] .
  * Chael  _et al._ (2021) A. Chael, M. D. Johnson, and A. Lupsasca, ApJ 918, 6 (2021), arXiv:2106.00683 [astro-ph.HE] .
  * Paugnat  _et al._ (2022) H. Paugnat, A. Lupsasca, F. H. Vincent, and M. Wielgus, A&A 668, A11 (2022), arXiv:2206.02781 [astro-ph.HE] .
  * Cunningham (1975) C. T. Cunningham, ApJ 202, 788 (1975).
  * Lindgren  _et al._ (2011) 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 .
  * Fuglstad  _et al._ (2015) G.-A. Fuglstad, F. Lindgren, D. Simpson, and H. Rue, Statistica Sinica 25, 115 (2015), arXiv:1304.6949 [stat.ME] .