Rotation in Event Horizon Telescope Movies

The Event Horizon Telescope (EHT) has produced images of M87* and Sagittarius A*, and will soon produce time sequences of images, or movies. In anticipation of this, we describe a technique to measure the rotation rate, or pattern speed $\Omega_p$, from movies using an autocorrelation technique. We validate the technique on Gaussian random field models with a known rotation rate and apply it to a library of synthetic images of Sgr A* based on general relativistic magnetohydrodynamics simulations. We predict that EHT movies will have $\Omega_p \approx 1^\circ$ per $GMc^{-3}$, which is of order $15\%$ of the Keplerian orbital frequency in the emitting region. We can plausibly attribute the slow rotation seen in our models to the pattern speed of inward-propagating spiral shocks. We also find that $\Omega_p$ depends strongly on inclination. Application of this technique will enable us to compare future EHT movies with the clockwise rotation of Sgr A* seen in near-infrared flares by GRAVITY. Pattern speed analysis of future EHT observations of M87* and Sgr A* may also provide novel constraints on black hole inclination and spin, as well as an independent measurement of black hole mass.


INTRODUCTION
The Event Horizon Telescope (EHT) has imaged the black hole Sagittarius A* (Sgr A*) at the heart of our own galaxy (Event Horizon Telescope Collaboration et al. 2022a) and the black hole M87* at the center of M87 (Event Horizon Telescope Collaboration et al. 2019a) at event horizon scale resolution. These images were made by combining data from an array of radio telescopes using a technique called very long baseline interferometry (VLBI). For M87*, key science results include a mass measurement that is consistent with estimates based on stellar kinematics (Gebhardt et al. 2011). For Sgr A*, key results include a mass measurement that is consistent with earlier, more-precise measurements based on individual stellar orbits (Schödel et al. 2002;Ghez et al. 2003Ghez et al. , 2008Do et al. 2019;GRAV-ITY Collaboration et al. 2019, 2020a. Interpretations of EHT data have relied heavily on time-dependent general relativistic magnetohydrodynamics (GRMHD) models, which are remarkably consistent with the data (Event Horizon Telescope Collaboration et al. 2019bCollaboration et al. , 2021Collaboration et al. , 2022bWong et al. 2022). In M87*, GRMHD models predicted (Event Horizon Telescope Collaboration et al. 2019b) that the angle between the brightness maximum on the ring and the large-scale jet in M87* observed in 2017, ∼ 150 • , was an outlier, and that an angle closer to ∼ 90 • would be more frequently observed. This is consistent with data from other epochs (Wielgus et al. 2020). In Sgr A*, however, GRMHD models predict a source-integrated variability that is a factor of two larger than observed (Event Hori- , and a plot of pressure averaged around the black hole's equatorial plane from the same time slice as the snapshot (right). In the right panel the region inside r = 2 GM c −2 is grayed out as it contributes relatively little emission. Notice that the nonaxisymmetric, time-variable structures form trailing spirals that are visible in both the snapshot and GRMHD pressure field. The field of view in the left and center panels is 20GM/(c 2 D), and the extent of the right panel is 20 GM c −2 .
zon Telescope Collaboration et al. 2022b; , focusing interest on the origins of variability in GRMHD models.
Variability is likely to become a focal point for EHT science. The EHT is developing the ability to revisit sources regularly, enabling movies of M87*, while also expanding its baseline coverage, enabling movies of Sgr A* (Doeleman et al. 2019;Johnson et al. 2019). What might movies reveal about both of the resolved EHT sources?
The hot spot model is a common starting point for understanding nonaxisymmetric variability. In the simplest version of this model, a hot spot moves freely on a circular orbit in the equatorial plane of the black hole (e.g. Broderick & Loeb 2006, Emami et al. 2023. Assuming emission arises near x ≡ Rc 2 /(GM ) ∼ 4, as it does in GRMHD-based models (see, e.g., Figure 4 of Event Horizon Telescope Collaboration et al. 2019b), then we expect the hot spot to orbit at the circular geodesic, or Keplerian, frequency Ω K = (GM/c −3 ) −1 (x 3/2 + a * ) −1 . For a face-on black hole with spin a * ≡ Jc/GM 2 = 0, Ω K ≈ 7(x/4) −3/2 degrees per GM c −3 . This frequency is an important point of comparison for variability in EHT movies.
GRMHD models do not show freely orbiting hot spots. Instead, they tend to show transient spiral features. Figure 1 shows a time-averaged image from a GRMHD model in the left panel next to a typical snapshot from the same model in the center panel. Evidently the nonaxisymmetric, time-dependent emission is concentrated in spiral features. The underlying plasma is subject to pressure gradient forces and magnetic forces, so the plasma need not move on geodesics. Strongly magnetized models (called magnetically arrested disks, or "MADs") tend to show rotation that is sub-Keplerian, while weakly magnetized models (called standard and normal evolution, or"SANEs") are closer to Keplerian. Radial velocities are typically close to the sound speed, particularly in models where the emission peaks inside the innermost stable circular orbit, in the so-called plunging region. The plasma motion is not well described by circular orbits.
The motion of the spiral features seen in Figure 1 may be detectable in EHT movies. In this paper, we define and evaluate the pattern speed Ω p , which is a measure of rotation in EHT movies. Our analysis is based on synthetic GRMHD data from the Illinois Sgr A* model library, which was run using KHARMA, an ideal nonradiative GRMHD code 1 and imaged with ipole (Mościbrodzka & Gammie 2018). The model library movies have an angular resolution of 0.5 µas and a time resolution of 5 GM c −3 between images. Each model comprises 3 × 10 3 images evenly spaced between time 1.5 × 10 4 to time 3 × 10 4 GM c −3 after their initialization with a magnetized torus. In Sgr A*, where GM c −3 20 s, the time between frames is 100 s and the total movie duration is 83 hr. For M87*, where GM c −3 9 hr, the time between frames is 2 days and the total movie duration is 15 yr. A more de- tailed description of how the library was made is provided in Event Horizon Telescope Collaboration et al. (2022b) and Wong et al. (2022).
This paper is organized as follows. Section 2 defines Ω p and introduces a methodology for measuring it in idealized synthetic image data. Section 3 applies this method to the Sgr A* model library and discusses the results. Section 4 provides a summary and describes next steps.

MEASURING PATTERN SPEED
The hot spot model discussed in Section 1 illustrates the difficulties in defining and measuring rotation in EHT movies. A single, equatorial, freely orbiting hot spot traces a complicated trajectory on the plane of the sky. Lensing can produce multiple images. Lensing has a particularly strong effect when the hot spot is seen edge on in the equatorial plane; then the brightest images trace a trajectory both above and below the black hole shadow. 2 Relativistic foreshortening and lensing make the apparent motion nonuniform; at modest inclination, the hot spot appears to move more quickly as it approaches the observer and more slowly as it recedes. Clearly there is a lot of potential information in EHT movies.
In this paper, we set aside this complexity and ask the most basic questions about the motion of brightness fluctuations on the ring. First, is it possible to determine if the fluctuations circulate clockwise or counter-clockwise on the sky? Second, is it possible to measure a characteristic rotation frequency, or pattern speed, Ω p ?
We begin by reducing the movie data to a manageable form. In each synthetic image, we sample the surface brightness T b on a circle defined by where the position angle (PA) parameterizes the location on the circle, i is the inclination angle between our line of sight and the angular momentum of the disk, and a * ≡ Jc/GM 2 is the dimensionless black hole spin. We use Bardeen's coordinates for x and y expressed in units of GM/(c 2 D) for a distance D to the source. This circle coincides with the critical curve (or shadow boundary) to first order in a * (Gralla & Lupsasca 2020). The synthetic images are smoothed using a Gaussian kernel with FWHM = 20 µas, the nominal EHT resolution. Since there are ∼ 3 resolution elements across Sgr A*'s ring, the brightness distribution sampled on the ring is insensitive to the precise radius and centering of the circle. Figure 2 shows the ring superposed over an example synthetic image.
Evaluating the surface brightness at each point on the ring over the entire duration of the movie yields T b = T b (PA, t), which we will call a cylinder plot. This cylinder plot is periodic in PA. The left panel of Figure 3 shows a cylinder plot for a fiducial model. Although we sample along a thin ring, blurring to the EHT's nominal resolution causes near-ring bright features to appear on the ring, giving the ring an effective thickness.
The cylinder plot shows characteristic diagonal bands. These thin bands appear near-vertical simply due to the aspect ratio. Each band corresponds to the movement of a bright feature around the ring. The features' orientation implies a net rotation toward positive PA (counterclockwise on the sky). The average slope of these features is the pattern speed Ω p , which we will measure using an autocorrelation analysis of the cylinder plot.

Normalization
The cylinder plot in Figure 3 is (1) brightest at PA ≈ 90 • , which corresponds to the approaching side of the accretion flow, and (2) exhibits fluctuations in source brightness over time, with a large dip in brightness near t 17, 500GM c −3 . An autocorrelation of the raw cylinder plot will be dominated by a few brightest features and will thus throw away information in low surface brightness features.
The bright feature in the cylinder plot near PA = 90 • is partially explained by Doppler boosting. 3 This brightness peak would appear even if the emission were axisymmetric. Both the time-averaged and fluctuating emission are amplified there. The asymmetry dominates the autocorrelation of the cylinder plot, downweighting information from low surface brightness PAs and reducing the accuracy of Ω p measurements. Assuming the signal-to-noise ratio is high, we would like to use the information available from fluctuations at all PAs.
The source brightness variations likewise amplify both the mean brightness and nonaxisymmetric fluctuations. Assuming that the signal-to-noise ratio is high, we would like to treat each snapshot on an equal footing.
3 Flow geometry and lensing also contribute to ring asymmetry.
In order to weight snapshots and PAs with different total fluxes more equally, we construct the cylinder plot using log(T b ). We then normalize by performing a mean subtraction along each time slice and each PA slice. The resulting cylinder plot has a mean of 0 along each column and row. This normalization procedure is independent of the order in which the mean subtraction is applied. The mean-subtracted logarithmic cylinder plot produces more accurate pattern speed measurements compared to a mean-subtracted linear plot (see Section 2.3 for accuracy tests). The right panel of Figure 3 shows the normalized cylinder plot, denotedT b (PA, t).

From Autocorrelation to Pattern Speed
Once the cylinder plot is normalized, we autocorrelatẽ T b (PA, t). Setting PA = PA + ∆PA and t = t + ∆t, the dimensionless autocorrelation function ξ is where denotes an average, σ 2 is the variance ofT b , and F is the Fourier transform. Notice thatT b is periodic in PA but not in t. We do not apply an explicit window function in time. The discontinuity that results from joining the beginning and ending of the time series with periodic boundary conditions at each PA produces power at high temporal frequencies that does not affect our analysis. Figure 4 shows the autocorrelation function for our fiducial model. Only the central part of the correlation function is shown. The tilt of the correlation function suggests a pattern speed Ω p ∼ 1 • per GM c −3 .  The pattern speed Ω p can be measured using second moments of ξ. In this subsection, for clarity, we use φ ≡ ∆PA and τ ≡ ∆t for the arguments of ξ. The relevant second moments are then The domain of integration will be specified below. We define Ω p by applying a shear transformation to the correlation function, and the integration region, until the off-diagonal moment vanishes. That is, we define φ = φ − Ω p τ , and adjust Ω p so that M φ τ = 0. Then = Ω p M τ τ .
The second equality follows from the definition of φ (the Jacobian of the shear transformation is 1; notice that the domain of integration must be transformed as well). The third equality follows from the definition of the moments. The final equality follows if Ω p is defined so that M φ τ = 0. Thus which is evidently dimensionally correct. The domain of integration should be set to maximize accuracy of the estimate of Ω p . A limited number of independent frames is used to estimate ξ, which introduces noise in ξ. The relative uncertainty in ξ increases away from the origin, and outside a few correlation lengths ξ is completely dominated by noise. If the domain of integration is too large then the moments are dominated by noise. Near the origin, however, pixelation of ξ also introduces errors. If the domain of integration is too small, then accuracy is lost. With these considerations in mind, we choose to integrate over a region with ξ > ξ crit . We set ξ crit = 0.8 to maximize measurement accuracy in our test problems (see Section 2.3) and minimize outliers in a survey of Ω p over the GRMHD model library.
To summarize, Ω p is estimated using the following procedure. Beginning with a high angular resolution synthetic movie: (1) smooth each frame to the nominal EHT resolution using a Gaussian kernel with FWHM = 20 µas; (2) sample the ring specified by Equation (1) in these frames to obtain the cylinder plot T b (PA, t); (3) take the log and mean subtract the cylinder plot to obtainT b (PA, t) (see Section 2.1); (4) calculate the correlation function ξ; (5) evaluate the moments of ξ at ξ > ξ crit ; and (6) calculate Ω p using Equation (11). 4

Verification
As a first test of the procedure, we created three movies containing a superposition of transient hot spots moving with constant angular frequency of either Ω hs = 1.23, 2.72, or 3.14 • per GM c −3 near the photon ring radius. The procedure recovers Ω p = Ω hs to within 4.3%.
As a second test of the procedure, we produce mock cylinder plots with a known pattern speed using sheared Gaussian random fields f . We begin with an unsheared random field f u (PA, t) with a Matern power spectrum Here m is the angular Fourier coordinate, ω is the temporal Fourier coordinate, and m o and τ o are constants. The power spectrum is sheared by setting P s = P u (m, ω + mΩ s ). Then a realization of the sheared field is generated in the Fourier domain from P s and transformed back to a realization in coordinate space f (PA, t), which is our mock cylinder plot. Calculating Ω p for 500 realizations of f with different shear rates Ω s , we are able to assess the accuracy of our measurement and the effects of pixelation. For mock cylinder plot parameters that are similar to those in the model library, we recover Ω p = Ω s with a root mean squared error of 2.8%.

PATTERN SPEEDS IN THE SGR A* MODEL LIBRARY
The Illinois component of the Sgr A* model library has four parameters: the magnetic flux (MAD or SANE), the inclination angle i, black hole spin a * , and the electron temperature parameter R high . In our convention, i is the angle between the line of sight and the accretion flow orbital angular momentum vector. Models with a * > 0 have prograde accretion flows, and models with a * < 0 have counterrotating, or retrograde, accretion flows. All models are assumed to have spin parallel or antiparallel to the accretion flow angular momentum; they are untilted. The electron temperature is set using the R high model, in which the ratio of the ionto-electron temperature varies smoothly from 1 where β 1 to R high where β 1. For more details on this prescription, see Wong et al. (2022), Equation (22). We have measured Ω p across the entire model library. Table  2 in the appendix lists Ω p for all models.

Sub-Keplerian Pattern Speeds
Our first main finding is that Ω p is small compared to what one would expect in a Keplerian hot spot model. The largest value measured in the entire library is 2.60 • per GM c −3 , and a more typical value is 1 • per GM c −3 (see Table 1). Thus Ω p Ω K (r = 4 GM c −2 )/7. This is small compared to what one would expect for hot spots orbiting freely close to the radius of peak emission.
The plasma is not orbiting freely, however, and instead exhibits pressure-driven and magnetic field-driven velocity fluctuations. In MAD models, the azimuthal fluid velocity is strongly sub-Keplerian. Could this explain the low pattern speeds? Figure 5 shows the orbital frequency of the underlying plasma u φ / u t , as seen by a distant observer in spherical Kerr-Schild coordinates. The figure shows the time-and longitude-averaged mean for both MADs and SANEs with a band indicating one standard deviation around the mean. The measured pattern speeds, which are shown as dashed lines spanning the principal emission region, are well below the fluid velocity. In SANE models, the azimuthal fluid velocity is indistinguishable from Keplerian. Apparently sub-Keplerian azimuthal fluid velocities do not provide a consistent explanation for the low pattern speeds.
The pattern speed seen in movies is, however, similar to the pattern speed measured in gas pressure (p g ) fluctuations. We measured p g (r = 3 GM c −2 , θ = π/2, φ, 1.5×10 4 < tc 3 /(GM ) < 3×10 4 ) to make a cylinder plot of gas pressure in the fiducial model. Figure 6 shows the raw cylinder plot and, after normalization, the correlation function. The pattern speed for the gas pressure fluctuations is Ω p = 1.1 • per GM c −3 , while the pattern speed for the synthetic images is Ω p = 1.1 • per GM c −3 . We find similar results in other models. It seems the sub-Keplerian pattern speeds are not an artifact of the low angular resolution in the images, radiative transfer, or lensing: they result from sub-Keplerian pattern speeds in the accretion flow itself.
The online animated Figure 7 shows the evolution of the image at both the full resolution of our synthetic images and the nominal EHT resolution. The clock hand in the animation moves at the pattern speed Ω p for this model, measured on the ring in the blurred images. The full resolution animation shows that the pattern speed is tracking narrow, trailing spiral features. The spiral features, like nearly all emission in MAD models, arise close to the midplane (see Figure 4 of Event Horizon Telescope Collaboration et al. 2019b).
The pattern speed associated with brightness fluctuations in the images and gas pressure fluctuations in the GRMHD models is slow compared to both the Keplerian speed and the azimuthal speed of the plasma. It must therefore be measuring a wave speed in the plasma. Since we see emission from only a narrow band in radius in these images (see the images in Figure 7 and the emission map in Figure 4 of Event Horizon Telescope Collaboration et al. 2019b), we must be seeing the azimuthal phase speed of the wave.
The underlying wave field is a combination of linear and nonlinear (shock) excitations. For simplicity, let . An unnormalized cylinder plot for the gas pressure at the midplane at r = 3 GM c −2 (left) in our fiducial GRMHD simulation (MAD, a * = 0.5), and its autocorrelation function ξ (right), where the blue line shows Ωp = 1.5 • per GM c −3 . This pattern speed is measured from the region with ξ > 0.3, which is inside the black contour. The vertical coordinate φ is the longitude in spherical Kerr-Schild coordinates.
us adopt a linear, hydrodynamic model, with a wave ∝ exp(ikr + imφ − iωt). In the tight winding (WKB) approximation, assuming that the disk is circular and Newtonian with orbital frequency Ω, the well-known dispersion relation is (ω −mΩ) 2 = Ω 2 +c 2 s k 2 . Here c s is the sound speed. The azimuthal component of the phase velocity is then ω/m = Ω ± (Ω 2 + c 2 s k 2 ) 1/2 /m. The phase velocity can thus be made small for the negative root and an appropriate choice of k and m. If the wave is trailing, as seen in the simulations, then the low phase velocity waves are ingoing. A nonlinear version of this argument is presented in Spruit (1987), which demonstrates the existence of stationary (zero pattern speed!) shocks in Keplerian disks.
Ingoing waves are plausibly excited by turbulence at larger radii. The pattern speed would correspond to the orbital frequency of the underlying plasma at the excitation radius (i.e. the corotation radius, where each dashed line of Figure 5 would intersect with the corresponding solid line). For a * = 0, the corotation radius ∼ 7 2/3 × 4 GM c −2 15 GM c −2 , well outside the region currently visible to the EHT. This suggests that disk fluctuations at radii outside the main emission region are important for determining variability.

Parameter Dependence
Our second main finding is that Ω p changes sign from i < 90 • to i > 90 • (recall that i is the angle between the line of sight and the accretion flow orbital angular momentum vector), so that the sign of Ω p signals the sense of rotation projected onto the sky. This stands in contrast to the time-averaged ring asymmetry, which follows the sign of a * (Event Horizon Telescope Collaboration et al. 2019b). It is therefore possible to distinguish between prograde and retrograde accretion in M87* by measuring both the ring asymmetry and Ω p (see Ricarte et al. 2022 for another technique for distinguishing retrograde accretion). This assumes the Sgr A* models considered here and M87* models have similar Ω p in degrees per GM c −3 , which is what we find in a sparse sampling of the M87* model library. Since the observed asymmetry shows that the spin vector in M87* is pointed away from Earth (Event Horizon Telescope Collaboration et al. 2019b), Ω p < 0 would imply prograde accretion and Ω p > 0 retrograde accretion.
Our third main set of findings concerns the dependence of Ω p on the model parameters. Fitting to Table 2 (in Section A), we find Ω p ≈ (1.5 + 0.4 a * − 0.2 ln R high ) cos i, SANE (12) in units of degrees per GM c −3 . The SANE fits have a root mean squared error of 0.31 • per GM c −3 , while the MAD fits have a root mean squared error of 0.14 • per GM c −3 . Systematic errors are discussed in Section 3.5. The worst-fitting models are generally SANE with a * > 0 and low R high . These models tend to have larger Ω p than otherwise expected. In Figure 8, we show the measured pattern speeds and the fits from Equations 13 and 12. The maximum SANE error is 1.15 • per GM c −3 larger than predicted, while the maximum MAD error is 0.42 • per GM c −3 larger. The fits in Equations (12) and (13) are consistent with the relatively slow rotation rate noted above (∼ 1 • per GM c −3 ) and with the sign of Ω p following the inclination (cos i dependence). Notice that for particles moving on a ring with angular frequency Ω, in flat space, the time-averaged apparent rotation rate would be Ω sgn(cos i).
The strongest dependences are on the inclination and mass (the mass dependence is in the units of Equations 12 and 13). In Sgr A*, the mass is known accurately from independent measurements (Schödel et al. 2002;Ghez et al. 2003Ghez et al. , 2008Do et al. 2019;GRAVITY Collaboration et al. 2019, 2020aEvent Horizon Telescope Collaboration et al. 2022c). A measurement of Ω p would thus constrain the inclination. In M87*, if we assume that inclination is determined by the large-scale jet, then | cos i| 1, so a measurement of Ω p would provide a distance-independent constraint on the mass.
The spin dependence of Ω p is weak but nonzero. It would seem to require accurate model predictions, lengthy observed movies, and careful interpretation of sparse interferometric data to use this dependence to constrain the spin.
We can estimate the uncertainty of these inclination, mass, and spin constraints using a probability distribution for Ω p at each set of parameter values obtained from kernel density estimation. This incorporates the uncertainties in measuring Ω p in movies of duration comparable to the expected observations (see Sections 2.3, 3.4, and 3.6). A single Ω p measurement could constrain inclination with a standard deviation of ∼ 20 • , with slightly smaller errors when the source is edge on and slightly larger errors when the source is face on. For M87*, if we assume the angle of the large-scale jet aligns with the inclination of the accretion flow, then an Ω p measurement would provide a distance-independent mass constraint with 1σ error of ∼ 33%.
The spin is unconstrained by a single Ω p measurement. Instead, the sign of Ω p will determine whether i > 90 • or i < 90 • . From there, the location of the peak brightness temperature will reveal whether the accretion is prograde or retrograde (see Figure 5 of Event Horizon Telescope Collaboration et al. 2019b). The spin is better constrained by making multiple Ω p measurements across different radii (see the right panel of Figure 9). These uncertainties do not account for the systematic errors in our models (see, e.g., Section 3.5), uncertainty in movie reconstructions from incomplete (u, v) coverage, or more informative priors.
Finally, notice that Ω p is nearly independent of R high . This suggests but does not prove that measurements of Ω p are insensitive to electron temperature assignment schemes. Ω p has a stronger R high dependence in SANE models, plausibly due to the stronger R high dependence of the emission region latitude in SANEs compared to MADs (see Figure 4 of Event Horizon Telescope Collaboration et al. 2019b). This limited R high dependence is also consistent with the close connection between the pattern speed in the gas pressure and the pattern speed in the images noted in Section 3.1.  . Four frames from an animation of the unblurred (top) and blurred (bottom) Sgr A* fiducial model, with a clock arm rotating at the measured pattern speed at the critical curve angular radius. The example frames are separated at a cadence of 10 GM c −3 . In the online animation, we show an unblurred (left) and blurred (right) movie of the Sgr A* fiducial model from 15, 000 to 20, 000 GM c −3 . By eye, we can see spiral shocks moving around the ring at the mean measured pattern speed, with some variation about the mean. The Keplerian velocity is over 5× the pattern speed, which is clearly too fast to fit the mean angular phase velocity of these spiral shocks. The EHT may observe at 345 GHz, providing higher resolution than existing data, which are taken at 230 GHz. Long-baseline space VLBI observations may enable even higher resolution. It is natural to ask whether the measured pattern speed changes with angular resolution. We can assess this by revisiting our analysis and smoothing the movies to different resolutions, recalling that up to now we have used Gaussian smoothing with FWHM = 20 µas, appropriate to the EHT's nominal resolution at 230 GHz. The left panel in Figure 9 shows the dependence of |Ω p | averaged across face-on Sgr A* models. It seems Ω p is only weakly dependent on resolution.

Dependence on Resolution and Sampling Radius
Our standard procedure for measuring Ω p samples brightness temperature fluctuations on a circle with angular radius b o = √ 27GM/(c 2 D), per Equation (1). This is close to the critical curve where the source is brightest and fluctuations are easiest to measure. What a * 0.94 a * 0.5 a * 0 a * -0.5 a * -0.94 Figure 9. The magnitude of Ωp averaged across all face-on models (i = 10 • , 170 • ) from time 15 × 10 3 to 20 × 10 3 GM c −3 , shown at different angular resolutions (left) and different radii (right). On the left, the nominal angular resolution corresponds to the FWHM of the Gaussian kernel used to smooth the synthetic images. The blue triangle corresponds to an angular resolution equivalent to 230 GHz observations, where the majority of this analysis was performed. The orange square and the green circle correspond to 345 and 690 GHz observations respectively, which exhibit small increases in Ωp and substantial increases in the variation over the model parameters, shown here as one standard deviation error bars measured over all MAD models. On the right, much of this analysis was done at the apparent critical curve angular radius bo; measurements at higher radii lead to a smaller average Ωp. The blue upward triangle, orange star, green circle, red box, and purple downward triangle correspond to a * = 0.94, 0.5, 0, −0.5, and −0.94 respectively.
would happen if fluctuations were measured at other impact parameters? The resulting "rotation curve" for nearly face-on models is shown in the right panel of Figure 9. When averaging over all inclinations, we find a fit consistent with Ω p 0.7(b o / √ 27) −1/2 , although the fit only covers a factor of 4 in b o , so the scaling is not strongly constrained. Evidently nonaxisymmetric fluctuations are not dominated by a single pattern speed at all impact parameters; instead there is spectrum of fluctuations with the dominant Ω p varying with b o . The impact parameter dependence also changes with spin, as seen in the right panel of Figure 9. Positive-spin (prograde) models show a greater change in Ω p with radius than negative-spin (retrograde) models.

Long-timescale Variability
We have checked for consistency of Ω p over time by subdividing our movies into three subintervals of duration 5 × 10 3 GM c −3 and measuring Ω p in each one. We find that the standard deviation between subintervals, averaged over all models, is ∼ 0.1 • per GM c −3 . Analysis of shorter-duration subintervals can be found in Section 3.6.
The variation between subintervals is larger for SANE models than MAD models, and for models with a * = 0.94. The variation is smaller for models with R high = 1. This long-timescale sample variance sets a fundamental limit on how accurately Ω p can be measured.

Light, Fast and Slow
Our model images are generated using the "fast light" approximation, which freezes the model on a single time slice and then ray traces through that time slice. Fast light is used because it is simple and the code is easy to parallelize. However, the fast light approximation fails to accurately represent changes in the source that occur on the light-crossing time. Short-timescale variations can be accurately captured by ray-tracing through an evolving GRMHD model, a procedure known as "slow light." Does the fast light approximation compromise our estimates of Ω p ? Figure 10 shows normalized cylinder plots for the fast light and slow light versions of a moderately inclined, high-spin model where fast light might be expected to have difficulty (MAD, a * = 0.94, i = 60 • , R high = 40). The cylinder plots differ in detail, especially on short timescales: we find the fast light model has Ω p = 0.31 • per GM c −3 , while the slow light model has Ω p = 0.37 • per GM c −3 , an increase of 0.06 • per GM c −3 or 19%. This increase is not large enough to change our conclusions, but it does imply additional uncertainty in Equations (12) and (13) that cannot be accurately evaluated without a slow light model library.

Pattern Speed Measurements in Observations
Measurements of Ω p may be affected by the observing cadence, duration, and by limited (u, v)  respectively. We have also measured Ω p in the example model from Section 2 (MAD, a * = 0.5, i = 30 • , R high = 160) at a cadence of 5, 10, and 20 GM c −3 , finding Ω p = 1.1, 1.1, and 0 • per GM c −3 respectively. For both models, Ω p decreases with the cadence, as the fastest features are often short lived. Ω p is nearly independent of the cadence below a threshold of 10 GM c −3 . Beyond this threshold, the autocorrelation peak begins to drop off more steeply, and pixelation effects limit our accuracy (see the end of Section 2.2). Lowering ξ crit from 0.8 to 0.4, we find Ω p = 0.94, a 25% decrease, for a 20 GM c −3 cadence.
Observing runs are likely to have shorter duration and lower cadence than our GRMHD movies. We examined the short-timescale variability of Ω p across the model library when measured at 10 GM c −3 cadence and 300 GM c −3 duration by subdividing the full 1.5 × 10 4 GM c −3 span of each model into subwindows. The results are shown in Figure 11, which shows the mean and standard deviation of Ω p over subwindows. The average standard deviation is 0.32 • per GM c −3 and the root mean squared variation is 0.06 • per GM c −3 .

CONCLUSION
This paper proposes a method to measure the rotation of brightness fluctuations in synthetic movies of EHT sources. We start by measuring surface brightness near the photon ring as a function of PA and time to produce a so-called cylinder plot. Rotation manifests in the cylinder plot as features that brighten, change PA, and then fade away. After normalizing the cylinder plot, we calculate its autocorrelation and use second moments of the autocorrelation to measure the apparent pattern speed Ω p .
We ran this procedure over the entire Illinois Sgr A* image library, which covers a broad range of plausible configurations for the source plasma, and in every case found the near-horizon pattern speed to be strongly sub-Keplerian, with a mean magnitude of |Ω p | ≈ 1 deg/GMc −3 , which is only approximately oneseventh of the expected Keplerian orbital velocity. This phenomenon can plausibly be attributed to the azimuthal phase velocity of ingoing spiral shocks excited outside the region that produces the bulk of the emission. Low pattern speeds are a fundamental prediction of GRMHD models.
We also found that Ω p depends strongly on inclination, with Ω p changing sign as the inclination crosses 90 • . The pattern speed scales with mass and is only weakly dependent on the spin. The expected pattern speeds are summarized by the fits in Equations (12) and (13).
In M87*, a pattern speed measurement would constrain the black hole mass, independent of distance, assuming the accretion flow inclination matches that of the large-scale jet. Since the sign of the pattern speed follows the accretion flow and the asymmetry of the observed ring follows the black hole spin (Event Horizon Telescope Collaboration et al. 2019b), a measurement of the pattern speed in M87* can distinguish between prograde and retrograde accretion.
In Sgr A*, where the black hole mass is known to high accuracy from stellar orbit measurements but the accretion flow inclination is not known, a pattern speed measurement would constrain the inclination. Wielgus  Figure 11. Ωp found with a longer 10 GM c −3 cadence and a shorter 300 GM c −3 duration, across all MAD models (left) and SANE models (right). The mean of each set of 300GM c −3 subwindows corresponds to the y-value, while the error bars denote the standard deviation. If there were no variability in Ωp[300 GM c −3 ] from that found in the full 15, 000 GM c −3 window, all points would fall along the dotted green line. et al. (2022) measured the linear polarization of millimeter emission immediately following an X-ray flare and found an evolution consistent with clockwise motion on the sky (see also Vos et al. 2022). The GRAVITY Collaboration has measured astrometric motion during infrared flares of Sgr A*, and this motion is also consistent with clockwise rotation on the sky and an approximately face-on inclination (GRAVITY Collaboration et al. 2018Collaboration et al. , 2020b. Both of these measurements lead us to expect that Ω p < 0 in Sgr A*. Notice, however, that GRAVITY measured a super-Keplerian rate of rotation, corresponding to a pattern speed much faster than what we find in this work. The accuracy of pattern speed measurements is limited by angular resolution, movie frame rate, movie duration, and the fast light approximation. We found that a broad range of sampling cadences around 5 GM c −3 produces similar pattern speeds. Subdividing our synthetic movies into shorter-duration movies produces similar but slightly varying pattern speeds. In a single test model computed with both the fast light approximation and slow light (no approximation), the pattern speeds differ by 19%. Finally, the analysis in this paper uses a set of models with similar initial and boundary conditions. It would be interesting to measure the pattern speed in alternative models (e.g. Ressler et al. 2020) since the pattern speed may be uniquely sensitive to conditions at radii outside the emission region.
The results presented here suggest that EHT will be able to measure pattern speeds, and that this measurement will provide valuable parameter constraints for M87* and Sgr A*. Techniques will need to be developed, however, that work directly with data in the (u, v) domain, and which are robust to gaps in (u, v) coverage and irregularly spaced data. We leave the problem of optimally extracting pattern speeds from EHT data for future work.
the Natural Sciences and Engineering Research Council of Canada through a Discovery Grant.
We are grateful to George Wong for providing the slow light models used in Section 3.5. We thank Steve Balbus, Alex Lupsasca, Maciek Wielgus, George Wong, and the anonymous referee for comments that greatly improved this paper.
The main text summarized Ω p measurements in the Sgr A* library using fits (see Equations 13 and 12). This appendix provides Ω p and Ω p,FIT for each model in the Illinois component of the Sgr A* model library. The units for Ω p and Ω p,FIT are degrees per GM c −3 .