The Variability of the Black Hole Image in M87 at the Dynamical Timescale

The black hole images obtained with the Event Horizon Telescope ( EHT ) are expected to be variable at the dynamical timescale near their horizons. For the black hole at the center of the M87 galaxy, this timescale ( 5 – 61 days ) is comparable to the 6 day extent of the 2017 EHT observations. Closure phases along baseline triangles are robust interferometric observables that are sensitive to the expected structural changes of the images but are free of station-based atmospheric and instrumental errors. We explored the day-to-day variability in closure-phase measurements on all six linearly independent nontrivial baseline triangles that can be formed from the 2017 observations. We showed that three triangles exhibit very low day-to-day variability, with a dispersion of ∼ 3 ° – 5 ° . The only triangles that exhibit substantially higher variability ( ∼ 90 ° – 180 ° ) are the ones with baselines that cross the visibility amplitude minima on the u – v plane, as expected from theoretical modeling. We used two sets of general relativistic magnetohydrodynamic simulations to explore the dependence of the predicted variability on various black hole and accretion-ﬂ ow parameters. We found that changing the magnetic ﬁ eld con ﬁ guration, electron temperature model, or black hole spin has a marginal effect on the model consistency with the observed level of variability. On the other hand, the most discriminating image characteristic of models is the fractional width of the bright ring of emission. Models that best reproduce the observed small level of variability are characterized by thin ring-like images with structures dominated by gravitational lensing effects and thus least affected by turbulence in the accreting plasmas.


INTRODUCTION
The Event Horizon Telescope (EHT) has produced horizon scale images of the black hole in the center of the M87 galaxy using Very Long Baseline Interferometry (VLBI) at a wavelength of 1.3 mm (Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f).Together with upcoming studies of the black hole in the center of the Milky Way, Sgr A*, these images have opened up numerous avenues to studying physics at event-horizon scales related to both accretion processes and gravity.The structure of the accretion disks and jets observed around the black holes will improve our understanding of the behaviour of plasma in these environments (Event Horizon Telescope Collaboration et al. 2019e,f).The measurement of the sizes and shapes of the black-hole shadows gives insights into deviations of the black hole spacetime from the Kerr metric (Psaltis 2019;Psaltis et al. 2020).
Accretion flows, which make up the horizon-scale environments of black holes, are expected to be highly variable by nature.The radiation observed from these regions at 1.3 mm is dominated by synchrotron emission from electrons which are heated and accelerated by magnetohydrodynamic (MHD) turbulence (Event Horizon Telescope Collaboration et al. 2019e).The accretion flow is also expected to evolve at the dynamical time scales near the innermost stable circular orbit (ISCO), which range from 4 to 56 minutes for Sgr A* (for M BH = 4.148 × 10 6 M ; see Gravity Collaboration et al. 2019) and 5 to 61 days for M87 (depending on the black hole spin; Bardeen et al. 1972).The timescales suggest that there * NSF Astronomy and Astrophysics Postdoctoral Fellow can be substantial variability in the source structure over a single observing night for Sgr A*, which makes it difficult to obtain a static image.However, we expect substantial source variability in M87 only for observations that span about a week or longer.Evidence for long-term variability in the image structure of M87 has been reported recently, based on VLBI observations spanning nearly a decade (Wielgus & The Event Horizon Telescope Collaboration 2020).
The EHT observed M87 in April 2017 for four nights with a maximum separation of six calendar days between observations.For a black hole of the inferred mass (M BH = 6.5 × 10 9 M ; Gebhardt et al. 2011a;Event Horizon Telescope Collaboration et al. 2019f), six days correspond to 16 GM/c 3 .This is approximately equal to ∼ 3.5 dynamical timescales at the ISCO for a maximally spinning black hole.Over this observation span, only a limited amount of variability in source structure was inferred between days, reflected as a small change in the thickness or azimuthal brightness distribution of the observed bright emission ring (Event Horizon Telescope Collaboration et al. 2019d).
Our goal in this paper is to use General Relativistic Magneto-Hydrodynamic (GRMHD) simulations of accretion flows, followed by General Relativistic ray tracing and radiative transfer, with parameters that are relevant to the M87 black hole, in order to explore the short-term variability in the source structure as it is imprinted in the interferometric observables.The EHT, being an interferometer, measures complex visibilities, which denote the two dimensional Fourier transform of the source image on the sky.The measurements are typically decomposed into visibility amplitudes and phases (Thompson et al. 2017).Image reconstruction techniques enable us to use these variables and create the map of the source brightness on the plane of the sky (see Event Horizon Telescope Collaboration et al. 2019d).
The visibility measurements in VLBI are affected by timedependent atmospheric and instrumental errors of the telescopes, particularly challenging for the short radio wavelengths (Thompson et al. 2017).The custom-built EHT data reduction pipeline addresses the specific requirements of the millimeter VLBI (Event Horizon Telescope Collaboration et al. 2019c).The pipeline uses Atacama Large Millimeter/submillimeter Array (ALMA) as an anchor station, leveraging its extreme sensitivity for the calibration of the entire EHT array (Blackburn et al. 2019,Janssen et al. 2019).As a result, the measured visibilities indicate sufficient phasestability to enable long coherent averaging and building up the high signal to noise ratio.Remaining station-based errors can be modelled as complex gains multiplying the visibilities.In order to eliminate the effects of these errors in the data, one can construct closure quantities with amplitudes and phases (Jennison 1958, Kulkarni 1989, Thompson et al. 2017, Blackburn et al. 2020).For a set of three baselines forming a closure triangle, the closure phase is defined as the sum of the measured visibility phases in each of these baselines.This cancels out the station-based phase measurement errors for each station.As a result, closure phases are optimal and robust probes of the intrinsic structure of the source.
Closure phases have indeed proven to be a powerful tool to study and quantify variability in synthetic data from GRMHD simulations.Medeiros et al. (2018) studied the dependence of variabililty in closure phases on the orientations and baseline lengths of the closure triangles (see also Roelofs et al. 2017).Triangles that involve large baselines (that probe small length scales in the source image) were shown to have a high degree of variability in closure phases.On the other hand, triangles that involve small baselines (that probe the overall size of the ring and source structure) exhibit less variability.The exception to this rule are triangles involving a baseline close to a deep visibility minimum.The visibility phases in regions of visibility minima are extremely sensitive to minor changes in source structure.This localization of variability in the Fourier space is vital in understanding variability in GRMHD simulations and in observations.
In this paper, we first present a data-driven model to quantify the variability of observed closure phases in the various linearly independent triangles of the M87 observations across the six days of observations in 2017.We apply this algorithm to the 2017 EHT data on M87 and identify three closure triangles that show a remarkably small degree of variability.We then explore the degree of variability produced in a large set of GRMHD simulations, with different magnetic field configurations and prescriptions of the plasma physics, and understand the effects of various model parameters on the variability in the models.Finally, we compare the predictions of the GRMHD simulations to the observations and discuss how the latter constrain the physical properties of the accretion flow in M87.

CLOSURE PHASE OBSERVATIONS OF M87
The EHT measures complex visibilities, given by where I(x, y) represents the total intensity of radiation at a given spatial location (x, y) on the image plane and (u, v) denote the Fourier frequencies corresponding to the (x, y) coordinates.However, due to atmospheric and instrumental effects, the measured visibilities do not represent the actual visibilities of the source.The measured visibilities are related to the source visibilities through complex "gains".The measured visibility between the i-th and the j-th telescopes, V i j , can be written as where g i and g j are the complex gains associated with the two telescopes and the star superscript indicates complex conjugates.We define the closure phase as the argument of the complex bispectrum along a closed triangle of three telescopes (Jennison 1958).The closure phase φ CP i jk for a closure triangle formed by the i-th, the j-th, and the k-th telescopes is hence given by where φ i j denotes the visibility phase between the i-th and the j-th telescopes.Equation 3 shows that the measured closure phase in a triangle is equal to the actual closure phase represented by the source.
The EHT observed M87 on four nights in 2017 -April 5, April 6, April 10, and April 11.The observations involved seven stations spread across five geographical locations (Event Horizon Telescope Collaboration et al. 2019b): the Atacama Large Millimeter Array (ALMA) and the Atacama Pathfinder Experiment (APEX) in Chile; the James Clerk Maxwell Telescope (JCMT) and the Submillimeter Array (SMA) in Hawai'i; the Large Millimeter Telescope (LMT) in Mexico; the Submillimeter Telescope Observatory (SMT) in Arizona; and the IRAM 30m telescope on Pico Veleta (PV) in Spain.On each night of observation, the rotation of the earth implies that each baseline traces out an elliptical path in the u − v space with the progression of time.
Figure 1 shows the u−v coverage for M87 during one of the observing days (left panel) as well as the visibility amplitudes measured as a function of baseline length on the same day (right panel).Two deep minima in visibility amplitude can be seen at baseline lengths of ∼ 3.4Gλ and ∼ 8.3Gλ encountered in three baselines (LMT-SMA, SMT-SMA, PV-SMA) along the east-west orientation (indicated in blue).In contrast, the minimum encountered by the ALMA-LMT baseline is only marginally deep (right panel of Figure 1).
We construct closure phase quantities using all the baselines shown in Figure 1.However, closure triangles with two telescopes at the same geographical location yield trivial closure phases, with deviations from zero arising due to changes in the large scale source structure alongside systematic and thermal errors (Event Horizon Telescope Collaboration et al. 2019c).While these triangles are useful to quantify systematic errors in closure phases, they do not capture any information about the source structure at horizon scales.Hence, we only choose triangles with telescopes at different geographical locations for our study.Since Chile and Hawai'i have two stations each, we pick the station with the higher signalto-noise ratio.We, therefore, consider five stations -ALMA, SMA, LMT, SMT, and PV, out of which six linearly independent closure phases can be constructed.We list the six triangles in Table 1.We choose the triangles such that all of them involve ALMA, which has the highest signal-to-noise ratio, and hence smaller uncertainties in closure phase measurements.
In Figure 2, we show the evolution of closure phases with time in Greenwich Mean Sidereal Time (GMST) coordinates over all four days of observations in three triangles that exhibit substantial day-to-day variability in closure phases over the span of the 2017 campaign (see also Event (ii) The ALMA-SMA-SMT triangle shows a swing of 180 • on each day of observation.In addition to this, the direction of this swing flips between the first two days and the later two days of observation.(iii) In the ALMA-PV-SMA triangle, it is difficult to quantify the degree of variability because of the reduced complex visibility amplitudes between PV and SMA baselines, and the low signal-to-noise ratio in some data points.However, a persistent structure in the evolution of closure phases is not observed.We note that each of these closure triangles involve one of the baselines (i.e., LMT-SMA, SMT-SMA, PV-SMA) that cross a deep visibility minimum in the E-W orientation.In addition, we note that in the ALMA-SMA-SMT triangle, the closure phases show little day-to-day variability until ∼ 19 GMST on each day.We attribute this behavior with the fact that the SMT-SMA baseline encounters the deep visibility minimum at ∼ 19 GMST.We plan to analyze the dependence of variability in the observations on the depth of the visibility minima in a future study.
The closure phases in the remaining triangles (i.e., ALMA-LMT-SMT, ALMA-PV-LMT, and ALMA-PV-SMT) show a persistent and continuous evolution with time during each day of observation (see Figure 3).This behaviour represents the evolution of closure phases as the various baselines trace their paths on the u − v plane following the rotation of the Earth.However, the same evolution with time is repeated in all days of observations; there is little scatter around the general trend set by the rotation of the baselines.The presence of substantial scatter from day to day would have been a signature of structure variability in the image.
In order to compare the observations to theoretical models, we quantify the degree of closure phase day-to-day variability observed in this last set of triangles.Upon examination of the closure phases in the maximal set of closure triangles, we find that outside the set of linearly independent triangles that we choose, only one triangle (LMT-SMT-PV) shows a low   7) in the low-variability triangles degree of day-to-day variability.We conclude that it is optimal and sufficient to quantify the day-to-day variability in the low-variability triangles that we choose in Table 1, since the LMT-SMT-PV triangle can be obtained by a linear combination of these three triangles.
We also note that the comparison between observations and models could, in principle, be carried out for both high variability and low variability triangles.However, as we will discuss in § 4, images from theoretical models uniformly show a high degree of phase variability in baselines that cross a visibility minimum and, as a result, a comparison with the high-variability triangles in the data turns out not to have much discriminating power between various models.The degree of variability exhibited on the low-variability triangles identified in the data, on the other hand, is a more significant challenge to the GRMHD models and proves to be a useful tool to distinguish between them.
The quantification of day-to-day variability in the set of three low-variable triangles is not straightforward, because closure phases evolve with time on each observing day but the individual scans on each day are not aligned at the exact same times.As a result, a difference in closure phase measured in two scans separated by almost (but not exactly) one sidereal day incorporates both the change due to the slightly different locations on the u − v plane probed by the two scans and the change due the structural changes in the image.
In order to disentangle the two sources of variability, we employ a data-driven analytic model for the evolution of closure phase with time during a single day of observation.This is meant to capture the change in closure phases introduced by the changing location on the u − v plane of the closure triangles.We then compare the overall change of the parameters of this model from day to day, in order to quantify the variability due to structural changes of the image.
For our data-driven model, φ model (t), we use a second order polynomial in the observation time t (in GMST coordinates), given by (4) In order to account for potential structural changes from day to day, we allow for an intrinsic Gaussian spread in the zeroth order coefficient of this polynomial (c 0 ), with a standard deviation of σ c0 .In this prescription, we do not assign any physical significance to the polynomial coefficients.Since the maximum separation between any two nights of observations is six days, σ c0 acts as a constraint on the overall change in closure phase at a given time in a given triangle across six days.
Using this model, we define the likelihood of making a particular closure phase measurement using a Gaussian mixture model as where φ i is the closure phase in a triangle at an observation time t i and σ i denotes the uncertainty in φ i .Note that we have assumed here Gaussian statistics for the measurements of closure phases.This assumption is not expected to introduce any significant changes in our quantitative results since the triangles we will be applying this method to are characterized by high signal-to-noise measurements.Applying this mixture model to data with large errors would require using an appropriate distribution for the closure-phase errors and their correlations (see Christian & Psaltis 2020).
The likelihood that n closure phase observations in a given triangle across all four days of observation follow the model φ model (t) is Assuming flat priors in the polynomial coefficients and in σ c0 , we map the parameter space using Markov-chain Monte Carlo (MCMC) sampling.
Figure 3 shows the posteriors over two of the model parameters, c 0 and σ c0 , for each of the three triangles that exhibit low day-to-day variability.The left panels show the 68% and 95% contours of the posteriors; the right panels show the data-driven model with the most likely value of the parameters, overplotted on the closure phase observations in these triangles.The most likely standard deviations of variability across the 6 days of observations are ∼ 3-5 • .This is comparable to the systematic error in closure phase observations of 2 • as reported in Event Horizon Telescope Collaboration et al. ( 2019c) and is remarkably small.In this paper, we analyze the GRMHD simulations that were run using harm (see Gammie et al. 2003 andNarayan et al. 2012).The simulations include SANE and MAD configurations of the magnetic fields and different black hole spins with co-rotating or counter-rotating accretion disks.The flow properties, as calculated in the GRMHD simulations, scale with the mass of the black hole and, therefore, the latter does not enter the calculation.The dependence of the electron temperature on the local value of the plasma β parameter aims to capture the effects of sub-grid electron physics that are not resolved in GRMHD simulations (see Chan et al. 2015).The accretion flows in the simulations are modeled as collisionless plasmas, where the ratio of ion temperature to electron temperature (R ≡ T i /T e ) is given by Here, β ≡ P gas /P mag is the ratio of gas pressure (P gas ) to the magnetic pressure (P mag ).This prescription aims to simulate the effects of sub-grid electron heating, using this simple, physics-motivated parametric form (Mościbrodzka et al. 2016, andEvent Horizon Telescope Collaboration et al. 2019e).
General relativistic ray tracing and radiative transfer is performed on the simulated system at the wavelength of 1.3 mm in order to generate the images of the black hole using the different snapshots of the simulations.For radiative transfer, the black hole mass (M BH ) and the overall electron number density scale in the accretion disk (n e ) set a scale for the mass and the optical depth in the accretion disk, which determines the mass accretion rate and the flux of emission in the source.It follows from the scaling properties of the emissivity of the accretion disk at 1.3 mm and from the natural scaling of the GRMHD simulations (see Chan et al. 2015) that n e and M BH are degenerate quantities with the resulting images and, therefore, the various interferometric observables, depending only on the product M BH n 2 e (see a derivation and discussion in Appendix A).
We use two sets of GRMHD models for studying closure phase variability: (i) The first set of models (henceforth Set A), is aimed at understanding the dependence of the closure-phase variability on the properties of the accretion flow, i.e., MAD vs. SANE magnetic field configuration, the plasma prescription, and the electron number density scale (n e ).In particular, they include SANE and MAD models for a single spin, three different values of R high , and five different values of n e (defined at a fiducial mass of M BH = 6.5×10 9 M ).For this value of the blackhole mass, this range allows us to explore compact fluxes that are up to a factor of a few above and below the nominal value of 0.5 Jy.Each model is used to generate 1024 images with a time cadence of 10 GM/c 3 , amounting to a total of ∼ 30000 images in this set.Ray tracing and radiative transfer calculations for this simulation set is carried out using GRay.
(ii) The second set of models (henceforth Set B) explore a comprehensive set of black hole spins and R high of the accretion plasma, which are used in Event Horizon Telescope Collaboration et al. (2019e) in the interpretation of the EHT observations.The electron number density scale (n e ) in this simulation set is tuned a priori such that the total flux of radiation at 1.3 mm is equal to 0.5 Jy for a fiducial mass of M BH = 6.2 × 10 9 M .Each model in this simulation set is used to generate 600-1000 images at a time cadence of 5 GM/c 3 , amounting to ∼ 50000 images in this set.Ray tracing and radiative transfer calculations for this simulation set is carried out using ipole.
The ray tracing and radiative transfer tracing calculations on the GRMHD simulations performed using ipole and GRay ignore the effects of finite light travel time.Bronzwaer et al. (2018) showed that the affect of the approximation on the lightcuve of flux density in a simulation is restricted to only a few percent.Depending upon the nature of spatial correlations in structural variability in the simulation, the computed images may carry either time-correlated or timeuncorrelated features.These features will have an effect on the amplitude of variability inferred.In this study, we ignore these effects in inferring the variability in the images.
The parameter space studied in sets A and B are listed in Table 2.Although the mass of the black hole M BH does not enter as an independent parameter in the calculation of the black-hole images, it affects our interpretation of the observations in the following ways: (i) The characteristic time unit (τ ) in the GRMHD simulations is set by the mass of the black hole as τ = GM BH /c 3 .A higher mass of the black hole implies that an observation span of 6 days would translate to a shorter elapsed time in units of τ .
(ii) The angular size of the image as seen by a distant observer depends on the black-hole mass (∆θ ∼ λ/M BH ), which implies that the EHT probes different regions in the Fourier space for black holes of different masses.A higher mass of the black hole translates to the EHT probing structures at smaller length scales.
The mass of M87 has been measured using two methods prior to the EHT observations (Event Horizon Telescope Collaboration et al. 2019e).Gebhardt et al. (2011b) measured a mass of M BH = 6.6 ± 0.4 × 10 9 M using stellar dynamics.Walsh et al. (2013) measured a mass of M BH = 3.5 +0.9 −0.7 × 10 9 M by studying emission lines from ionized gas.The GRMHD models in simulation sets A and B give images for which the ring sizes and widths correspond to masses in the range inferred by the stellar dynamics measurement of Gebhardt et al. (2011b) (see Event Horizon Telescope Collaboration et al. 2019f).Because of this, we allow the black hole mass to vary in the range 5 × 10 9 M ≤ M BH ≤ 9.5 × 10 9 M in both sets of the simulations.
The orientation of the black-hole spin axis has been constrained by long-wavelength observations of its jet to be at a position angle of ∼ 288 • North of East and at an inclination of 17 • with respect to the line of sight (see discussion in Event Horizon Telescope Collaboration et al. 2019f).We trace the Set A models at an inclination of 17 • .In Set B, we trace the positive spin models at an inclination of 163 • and the negative spin models at 17 • .
The simulations typically run for lengths that correspond to ∼ 10 4 GM/c 3 .Here, we only use images obtained from the relaxed turbulent state with a slowly varying mass accretion rate.This corresponds to simulation times in the range of 5 × 10 3 ≤ t ≤ 10 4 GM/c 3 .

CLOSURE PHASE VARIABILITY IN GRMHD SIMULATIONS
In order to analyze visibility phase and closure phase variability in GRMHD models, we first transform the GRMHD snapshots into the visibility space by performing a twodimensional Fourier transform on each of the images.We  c The electron number density is defined at a fiducial mass of 6.5 × 10 9 M for set A and 6.2 × 10 9 M for set B.
then compute the amplitude and phase from the complex visibilities at each location in the u − v plane.
Since the visibility phase is a directional quantity, we need to employ directional statistics to infer the standard deviation in a time-series.For a given time-series of n phases θ i , we compute its directional dispersion as where θ represents the circular mean, defined as In the limit of small deviations in the time-series from the circular mean θ, the directional dispersion can be approximated as where σ is inferred as a standard deviation of the time-series θ i for small fluctuations about θ.
We use this quantity to construct heat maps of phase variability in order to understand the dependence of variability on the u − v coordinates (see also Medeiros et al. 2017).We present an example of such a heat map for one of the simulations in set B in Figure 4.The left panel shows the mean visibility amplitude on the u − v plane computed across the simulation run and reveals deep visibility minima that are aligned with the spin axis of the black hole.The right panel shows a heat map of visibility phase dispersion.There are "hot regions" of visibility phase variability, i.e., localized regions in u − v space that show large dispersions.These regions lie primarily along the spin axis of the black hole and at baseline lengths that correspond to the deep minima in visibility amplitudes.In fact, a natural anticorrelation between phase variability and mean visibility amplitude is observed in all the GRMHD models, which was also explored in Medeiros et al. (2017) (see their Figures 3 and 4).
Based on this understanding, one expects the closure phase variability computed along various triangles to depend on the locations of the three vertices in the u − v space.In particular, the localized nature of variability implies that closure triangles that have baselines crossing the hot regions are expected to exhibit high variability in closure phases, whereas triangles with vertices in quiet parts of the u − v space should show a low level of variability.In Figure 4, we show examples of two such triangles: the ALMA-SMA-LMT triangle has a baseline on a hot region and is expected to be highly variable, while a low-variability closure triangle (ALMA-LMT-SMT) avoids the hot regions.
The azimuthal extent of the hot regions of phase dispersion in the u − v plane depends on the degree of symmetry of the underlying image.This is shown in Figure 7, where each column corresponds to a simulation from set A with a different electron density scale n e , and all other black hole and plasma parameters held fixed.Qualitatively, the two effects of increasing the electron density scale are a higher level of symmetry in the bright emission ring and an increase in its fractional width.The more symmetric images lead to a larger azimuthal extent in the minima in the visibility amplitudes (middle panels) and in the hot regions of phase dispersion (lower panels).The increased fractional width of the images lead to the locations of both the visibility amplitude minima and the hot regions in phase dispersion to appear at smaller baseline lengths.

Closure Phase Variability and Compliance Fractions
The phase dispersions discussed in the previous section were calculated across the entire span of the simulations, which corresponds to many dynamical timescales.In order to compare the simulations to the observed variability in M87, however, we need to calculate the degree of variability for a timespan of ∆t = 6 days, which is the longest separation between the four observations in 2017.Furthermore, because observations only yield closure phases, hereafter we focus on this quantity to facilitate comparisons.
We define the variance in closure phase φ cp (t) for a timespan of ∆t = 6 days using Equation (11) as In order to reduce the discretization noise when calculating this quantity, we perform a linear interpolation between simulation time-steps and construct a continuous function φ cp (t) representing the evolution of closure phases.We use fixed u − v locations taken at a median time over the night of observation, in order to construct this function.This ensures that we separate the evolution of closure phases due to the changing locations of the baselines from the structural variability of the images that we are interested in inferring.
In Figure 5, we show the evolution of closure phase with time during one segment of a simulation chosen from set B, that has a MAD configuration with parameters a = +0.94 and R high = 20.As is evident from this evolution, there are time spans where the closure phase shows little variability (e.g., the times between the green vertical lines) and others where there are substantial swings of order π over a short period of time that is comparable to the length of the observations.Given the low degree of variability observed for some trian-gles in M87 data (see § 2), this motivates us to define a metric that quantifies how common such periods are in a given model.
Figure 6 shows the standard deviation in closure phase variability computed from all 6-day segments in the simulation shown in Figure 5, for the three triangles that have observationally shown a small degree of variability.The dashed lines correspond to thresholds σ 1 , σ 2 , σ 3 of this observed variability.In Figure 5, we choose as the most likely values of σ c0 in the corresponding triangles as the thresholds (see Table 1).The green points indicate the occurrences of 6-day segments that are consistent with observations.We define the quantity F (σ 1 , σ 2 , σ 3 ) of a model as the fraction of 6-day segments that show a level of variability consistent with the observed one; i.e., the ratio of green points to the total number of points in this figure .Formally, F (σ 1 , σ 2 , σ 3 ) can be written as where P(σ sys ) is the prior in the systematic uncertainty σ sys and P model (σ 1 , σ 2 , σ 3 |σ sys ) is the three-dimensional distribution of circular standard deviations σ 1 , σ 2 , σ 3 computed in the three low-variability triangles using Equation 12.The observed variability, as quantified in § 2, is a combination of the true intrinsic variability of the source and the systematic errors in measurement.Because of our limited knowledge of the systematic uncertainties and their priors, we assume a form such that equation ( 13) can be written as i.e., simply assuming that the combined intrinsic degree of variability and the systematic uncertainties cannot exceed the inferred values.
We compute this compliance fraction for black hole masses in the range 5 × 10 9 M ≤ M BH ≤ 9.5 × 10 9 M and use the maximum number as the compliance fraction for a given GRMHD model.

COMPARISON OF GRMHD SIMULATIONS TO M87 OBSERVATIONS
In this section, we use the two sets of simulations we discuss in §3 to explore the dependence of the compliance fraction on the magnetic field configuration, electron temperature prescription in the plasma, and electron-density scale (set A), as well as on the black-hole spin (set B).These two sets are meant to serve as two slices in the large parameter space of simulations.
Table 3 shows the compliance fractions calculated for the 30 simulations in Set A. Overall, these compliance fractions are rather small, indicating that a lot of simulations show more frequent periods of high variability, even in the triangles that observationally do not exhibit large changes in closure phase across the 6 days of the campaign.The dominant distinguishing characteristic of the simulations that show a higher compliance fraction is the small value for the electron density scale.There is smaller difference in compliance fractions between the MAD and SANE simulations.We observe that for small values of R high , the dependence of compliance fraction on the fractional width is stronger.
We explored how the electron density scale affects the image structure and variability.To this end, we calculated the characteristic properties of the average image in each simulation, using the image-domain techniques discussed in Event Horizon Telescope Collaboration et al. (2019d) andEvent Horizon Telescope Collaboration et al. (2019f).We measured the fractional width of each average image δw f as the ratio of the width δw to the ring diameter d, i.e. δw f ≡ δw/d.We found the fractional widths of the images correlate strongly with the electron density scale, as shown in Figure 8.This is expected, given that increasing the electron density scale also increases the emissivity and optical depth in the accretion flow, as shown in the top panels of Fig. 7 (Chan et al. 2015).
The right panel of Figure 8 shows indeed that the compliance fraction is anticorrelated with the fractional width of the average images of the simulations.When the bright emission ring in an image is narrow, it traces the outline of the black hole shadow and, therefore, its shape and brightness distribution is determined primarily by gravitational lensing effects.On the other hand, when the bright emission ring is broad, variability due to the turbulence in the accretion flow introduces more apparent effects.The low level of variability In the simulations of Set B, the electron number-density scale was tuned to generate a fixed source flux of 0.5 Jy.Different configurations were explored with both rotating and counterrotating black-hole spins of various magnitudes.The compliance fractions are shown in Table 4 and are, overall, lower than those of the Set A simulations.As in the previous case, the MAD simulations have on average a higher compliance fraction than the SANE simulations and there is only a marginal dependence on the electron temperature parameter R high .However, there is a significant dependence on the spin of the black hole, with co-rotating, high-spin simulations producing the highest compliance fractions.
The table also shows constraints imposed on the models of Set B by other considerations, not originating from the EHT observations (Event Horizon Telescope Collaboration et al. 2019e).They include constraints related to (a) the simulations having achieved radiative equilibrium, (b) the predicted X-ray flux to be L X = (4.4± 0.1) × 10 40 erg s −1 , as measured during the nearly simultaneous observations with the Chandra X-Ray Observatory and the Nuclear Spectroscopic Telescope Array (NuSTAR), and (c) the jet power being in range 10 42 to 10 45 erg s −1 .Simulations with fast black-hole spins, which are consistent with these constraints, are also characterized by larger compliance fractions, giving preference to these models.

DISCUSSION
In this paper, we studied the variability of the black-hole image structure of M87 at timescales comparable to the fastest dynamical timescale near the horizon, as it is imprinted on the closure phases measured during the EHT 2017 observations.We identified three linearly independent closure triangles that exhibit a persistent evolution pattern of closure phases over the course of each night of observing but show little variation around this pattern across the six days of observations, namely, ALMA-LMT-SMT, ALMA-PV-SMT, and ALMA-PV-LMT.In other words, the closure phase evolution in these triangles follows set tracks determined by the rotation of the baselines on the u − v plane but shows very little scatter around these tracks from day to day.We quantified the degree of variability in these triangles to be ∼ 3 − 5 • .This inferred level of variability is comparable to the ∼ 2 • systematic error in measurement of closure phases (Event Horizon Telescope Collaboration et al. 2019c).
We also found three triangles that exhibit a high level of day-to-day variability, namely, ALMA-SMA-LMT, ALMA-SMA-SMT, and ALMA-PV-SMA.The change in closure phase across 6 days in a given location on the u − v plane for these triangles can be ∼ 90 − 180 • .We identified them as triangles with at least one baseline that encounters a deep visibility minimum on the u − v plane.This is in agreement with expectations based on theoretical models that reveal the presence of highly variable but highly localized regions on the u − v plane associated with the locations of these minima.
We used GRMHD simulations to explore the dependence of closure-phase variability on different model parameters and at different locations in the u − v plane.We found that the   most discriminating image characteristic of models in terms of the degree of closure phase variability is the fractional width of the ring of high intensity on the image.Models that best reproduce the observed small level of variability are those with thin ring-like images with structures dominated by gravitational lensing effects and thus least affected by turbulence in the accreting plasmas.
Among the models we explored, there is marginal difference between SANE and MAD simulations, which explore different magnetic field configurations in the accretion flows, with some preference for the MAD models.There is also a small dependence on black-hole spin, with high-spin corotating models showing the lowest level of day-to-day variability, in agreement with the observations.These findings demonstrate that the method we introduced to quantify day-to-day variability in closure phase data and compare it to models is a useful tool in exploring the origin of variability in horizon-scale images of black holes and in discriminating between models.
In order to calculate the various simulated images, we used the plasma properties from long GRMHD simulations.Nonradiative GRMHD simulations are invariant to a rescaling of the density by a factor M, as long as the magnetic field also is also rescaled by a factor of M 1/2 and the internal energy by a factor of M. In other words, B ∼ n 1/2 e and the Alfvén speed B/(m p n e ) 1/2 is the quantity that remains invariant under rescaling.Combined with the approximate expressions (A12)-(A13) derived above, this implies that the synchrotron emissivity and opacity evaluated using the simulations outputs are invariant to rescaling, as long as the product n e B α ∼ n 1+α/2 e remains constant.Finally, the integration of the transfer equation is performed on a coordinate system in which the distances are expressed in terms of the length scale set by the mass of the black hole, i.e., We drop the ν in subscript to indicate quantities calculated at ν = 230 GHz.From this last expression it is apparent that the electron number density n e and the black hole mass M BH are degenerate quantities in the solution of the radiative transfer problem at wavelengths where the dominant source of opacity is due to synchrotron processes, with a degeneracy in the product n 1+α/2 e M BH and with α 2.

Figure 1 .
Figure 1.(Left) : The u − v coverage traced by the EHT for the M87 observations on 2017, April 11, with the dashed circles indicating locations of the two observed visibility minima.The colors indicate different orientations of the u − v coverage.(Right): The visibility amplitude as a function of baseline length on the same day.
Horizon Telescope Collaboration et al. 2019c Section 7.3.2):(i) The ALMA-SMA-LMT triangle shows a change of 30 • − 60 • between the first two days and the later two days of observation.
baselines indicate the triangles which cross one of the deep visibility minima in the E-W orientation at 3.4 Gλ and 8.3 Gλ. b Baselines labeled 1, 2, and 3 in a triangle labelled as A-B-C correspond to the baselines A-B, B-C. and C-A respectively.c An estimated range of variation in the high-variability triangles, and the maximum likelihood value of the inferred Gaussian variability (Equation

Figure 2 .
Figure 2. The evolution of closure phase plotted with time for all four days of observation for closure triangles in which baselines are known to encounter regions of deep visibility minima.All the three triangles exhibit a high level of variability in closure phases across six days of observation.
3. GRMHD SIMULATIONS AND LIBRARYA large suite of GRMHD simulations has been generated to model the accretion flow around M87.The simulations have been performed using the algorithms harm (seeGammie et al. 2003), BHAC (Olivares Sánchez et al. 2018), KORAL (S ądowski et al. 2014), etc.The simulations have been initialized with two magnetic field configurations that led to different field structures in the flows: SANE (Standard and Normal Evolution, see Igumenshchev et al. 2003) and MAD (Magnetically Arrested Disk, see Narayan et al. 2012 and S ądowski et al. 2013).A comprehensive comparison of the GRMHD algorithms is presented in Porth et al. (2019).The set of observables that result from these simulations have been obtained by general relativistic ray tracing and radiative transfer algorithms, such as GRay (Chan et al. 2013) and ipole (Mościbrodzka & Gammie 2018).

Figure 3 .
Figure 3.The parameters and functional forms of the data-driven closure phase evolution model, applied to the closure triangles that exhibit low variability across the 6 days of 2017 EHT obsevations (see eq. 4).(Left): The posterior projected on the σc 0 − c0 parameter plane, where σ0 is a measure of the variability across the 6 days.The red dots indicate the most likely values.The contours indicate the levels containing 68% and 95% of the posterior probability.The dashed horizontal lines indicate the systematic error of 2 • in closure phase observations.(Right): Closure phase observations plotted as a function of time for the same triangles.The cyan band indicates the most likely value of the width of the model σc 0 .
u − v coordinates of the baselines are rotated to align the spin axis of the black hole at a position angle of 288 • East of North.b The positive spin models are ray traced at an inclination of 163 • and negative spin models at 17 • in set A. Set B models are ray traced at an inclination of 17 • .

Figure 4 .
Figure 4. Normalised mean visibility amplitude map on a logarithmic scale (left) and directional dispersion map (see Equation 9) of visibility phases (right) for the SANE, a = +0.5,Rhigh = 1 model from simulation set B. The black hole spin points upwards as indicated by the white arrow labeled as Ω.The EHT baseline coverage (on April 11) are shown by the white dots, and are rotated such that the spin points to a position angle of 288 • east of north (with north indicated by the white arrow labeled as N and the urot − vrot plane representing the rotated u − v plane).The cyan and the red triangles indicate the ALMA-LMT-SMT (at UTC 3:32:5.0003h / GMST 16:26:37.1531h) and ALMA-SMA-LMT (at UTC 5:00:5.0001h / GMST 17:54:51.6089h) triangles respectively.The second triangle includes the SMA-LMT baseline that crosses a visibility minimum, which translates to a high level of variability in closure phase introduced by structural changes in the image.The mass of the black hole is set to MBH = 7.5 × 10 9 M

Figure 5 .
Figure 5. Closure phases for three triangles plotted as a function of time for a particular segment of the same simulation (set B -MAD a = +0.94,Rhigh = 20).The black parallel lines indicate the time scale of 6 days and a standard deviation of ∼ 4 • .The section within the green dashed vertical lines indicates a region of low-variability, and hence counts towards the compliance fraction of the model.

Figure 6 .
Figure 6.Two dimensional scatter plots of inferred standard deviation in closure phases in three triangles across one simulation (simulation set B -MAD a = +0.94,Rhigh = 20).The horizontal and vertical dashed lines indicate the observational bounds for the standard deviation for each triangle, as described in § 2.across the 6 days observed in M87 argue in favor of images characterized by thin (∼ 20%) rings of emission.In the simulations of Set B, the electron number-density scale was tuned to generate a fixed source flux of 0.5 Jy.Different configurations were explored with both rotating and counterrotating black-hole spins of various magnitudes.The compliance fractions are shown in Table4and are, overall, lower than those of the Set A simulations.As in the previous case, the MAD simulations have on average a higher compliance fraction than the SANE simulations and there is only a marginal dependence on the electron temperature parameter R high .However, there is a significant dependence on the spin of the black hole, with co-rotating, high-spin simulations producing the highest compliance fractions.The table also shows constraints imposed on the models of Set B by other considerations, not originating from the EHT observations (Event Horizon Telescope Collaboration et al. 2019e).They include constraints related to (a) the simulations having achieved radiative equilibrium, (b) the predicted X-ray flux to be L X = (4.4± 0.1) × 10 40 erg s −1 , as measured during the nearly simultaneous observations with the Chandra X-Ray Observatory and the Nuclear Spectroscopic Telescope Array (NuSTAR), and (c) the jet power being in range 10 42 to 10 45 erg s −1 .Simulations with fast black-hole spins, which are consistent with these constraints, are also characterized by larger compliance fractions, giving preference to these models.

Figure 7 .
Figure 7.The average images (top), normalized average visibility amplitudes (middle), and directional phase dispersion (bottom) of the MAD, a = 0.9, R high = 20 simulation.Models with three electron number densities -ne = 1 × 10 5 cm −3 (left), ne = 5 × 10 5 cm −3 (middle), and ne = 1 × 10 6 cm −3 (right), are mapped in this figure.The white arrows labeled as Ω indicate the spin vector of the black hole, and the ones labeled as N indicate the north direction.The coordinates (urot, vrot) represent the u − v coordinates rotated such that the spin points to a position angle of 288 • east of north.It is evident that a model with lower ne exhibits emission dominant from the photon ring, and has a lower variability in visibility phases outside the first visibility minimum.The three triangles plane indicate the low-variability triangles -ALMA-LMT-SMT (cyan), ALMA-PV-SMT (green), and ALMA-PV-LMT (white) at UTC 3:32:5.0003h/ GMST 16:26:37.1531h.The colormap in visibility amplitude is in logarithmic scale with the maximum value is normalised to unity.

Figure 8 .
Figure 8. Left -The fractional width of the ring of emission in the average image of each simulation in set A. Right -The anticorrelation between fractional width and compliance fraction in the simulation set A.

Table 1 .
Closure Triangles used with Baseline Lengths and Inferred Variability in Closure Phases.

Table 3 .
Compliance fraction for Set A
c Satisfies Jet Power constraint.