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

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 (EHT Collaboration et al. 2019a, 2019b,2019c, 2019d, 2019e, 2019f). 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 behavior of plasma in these environments (EHT Collaboration et al. 2019e, 2019f). The measurements of the sizes and shapes of the black hole shadows give 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 that are heated and accelerated by magnetohydrodynamic turbulence (EHT Collaboration et al. 2019e). The accretion flow is also expected to evolve at the dynamical timescales near the innermost stable circular orbit (ISCO), which range from 4 to 56 minutes for Sgr A* (for M_BH = 4.148 × 10⁶ 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 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 EHT Collaboration 2020).

The EHT observed M87 in 2017 April 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⁹ M_⊙; Gebhardt et al. 2011; EHT Collaboration et al. 2019f), the six days correspond to 16 GM/c³. 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 (EHT Collaboration et al. 2019d).

Our goal in this paper is to use general relativistic magnetohydrodynamic (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 EHT Collaboration et al. 2019d).

The visibility measurements in VLBI are affected by time-dependent 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 (EHT 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 phase-stability to enable long coherent averaging and building up the high signal-to-noise ratio. Remaining station-based errors can be modeled 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, the 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 variability 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.

The EHT measures complex visibilities, given by

$\begin{eqnarray}&&V(u,v)=\int \int {dx}\ {dy}\ I(x,y)\ {e}^{-2\pi i({ux}+{vy})},\end{eqnarray} \tag{ 1 }$

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 ith and the jth telescopes, ${{\mathscr{V}}}_{{ij}}$, can be written as

$\begin{eqnarray}&&{{\mathscr{V}}}_{{ij}}={g}_{i}{g}_{j}^{* }{V}_{{ij}},\end{eqnarray} \tag{ 2 }$

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 ${\phi }_{{ijk}}^{\mathrm{CP}}$ for a closure triangle formed by the ith, the jth, and the kth telescopes is hence given by

$\begin{eqnarray}\begin{array}{rcl}{\phi }_{{ijk}}^{\mathrm{CP}} & = & \mathrm{Arg}[{{\mathscr{V}}}_{{ij}}{{\mathscr{V}}}_{{jk}}{{\mathscr{V}}}_{{ki}}]\\ & = & \mathrm{Arg}[{g}_{i}{g}_{j}^{* }{V}_{{ij}}{g}_{j}{g}_{k}^{* }{V}_{{jk}}{g}_{k}{g}_{i}^{* }{V}_{{ki}}]\\ & = & \mathrm{Arg}[| {g}_{i}{| }^{2}| {g}_{j}{| }^{2}| {g}_{k}{| }^{2}{V}_{{ij}}{V}_{{jk}}{V}_{{ki}}]\\ & = & \mathrm{Arg}[{V}_{{ij}}{V}_{{jk}}{V}_{{ki}}]\\ & = & {\phi }_{{ij}}+{\phi }_{{jk}}+{\phi }_{{ki}},\end{array}\end{eqnarray} \tag{ 3 }$

where ϕ_ij denotes the visibility phase between the ith and the jth 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, 6, 10, and 11. The observations involved seven stations spread across five geographical locations (EHT 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 Hawaii; the Large Millimeter Telescope (LMT) in Mexico; the Submillimeter Telescope Observatory (SMT) in Arizona; and the IRAM 30 m telescope in 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.4 and ∼8.3 Gλ encountered in three baselines (LMT-SMA, SMT-SMA, and 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).

Zoom In Zoom Out Reset image size

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 (EHT 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 Hawaii have two stations each, we pick the station with the higher signal-to-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 EHT Collaboration et al. 2019c, Section 7.3.2):

Zoom In Zoom Out Reset image size

(i) The ALMA-SMA-LMT triangle shows a change of 30°–60° between the first two days and the later two days of observation.

(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 involves one of the baselines (i.e., LMT-SMA, SMT-SMA, and 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 to 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 behavior 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.

Zoom In Zoom Out Reset image size

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 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 Section 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 phases 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 to 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 the 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

$\begin{eqnarray}&&{\phi }_{\mathrm{model}}(t)={c}_{0}+{c}_{1}t+{c}_{2}{t}^{2}.\end{eqnarray} \tag{ 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₀), with a standard deviation of ${\sigma }_{{c}_{0}}$. 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, ${\sigma }_{{c}_{0}}$ 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

$\begin{eqnarray}\begin{array}{rcl}{L}_{i}({c}_{0},{c}_{1},{c}_{2},{\sigma }_{{c}_{0}}) & \equiv & \left\{\displaystyle \frac{1}{\sqrt{2\pi ({\sigma }_{i}^{2}+{\sigma }_{{c}_{0}}^{2})}}\right\}\\ & & \times \,\exp \left[-\displaystyle \frac{1}{2}\displaystyle \frac{{\left\{{\phi }_{i}-{\phi }_{\mathrm{model}}({t}_{i})\right\}}^{2}}{\left({\sigma }_{i}^{2}+{\sigma }_{{c}_{0}}^{2}\right)}\right],\end{array}\end{eqnarray} \tag{ 5 }$

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

$\begin{eqnarray}&&L({c}_{0},{c}_{1},{c}_{2},{\sigma }_{{c}_{0}})=\prod _{i=1}^{n}{L}_{i}({c}_{0},{c}_{1},{c}_{2},{\sigma }_{{c}_{0}})\end{eqnarray} \tag{ 6 }$

$\begin{eqnarray}&&=\prod _{i=1}^{n}\left\{\displaystyle \frac{1}{\sqrt{2\pi ({\sigma }_{i}^{2}+{\sigma }_{{c}_{0}}^{2})}}\right\}\exp \left[-\sum _{i=1}^{n}\displaystyle \frac{1}{2}\displaystyle \frac{{\left\{{\phi }_{i}-{\phi }_{\mathrm{model}}({t}_{i})\right\}}^{2}}{({\sigma }_{i}^{2}+{\sigma }_{{c}_{0}}^{2})}\right].\end{eqnarray} \tag{ 7 }$

Assuming flat priors in the polynomial coefficients and in ${\sigma }_{{c}_{0}}$, we map the parameter space using the Markov-chain Monte Carlo sampling.

Figure 3 shows the posteriors over two of the model parameters, c₀ and ${\sigma }_{{c}_{0}}$, for each of the three triangles that exhibits 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 six days of observations are ∼3°–5°. This is comparable to the systematic error in closure-phase observations of 2° as reported in EHT Collaboration et al. (2019c) and is remarkably small.

A large suite of GRMHD simulations has been generated to model the accretion flow around M87. The simulations have been performed using the algorithms harm (see Gammie 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: Standard and Normal Evolution (SANE; see Igumenshchev et al. 2003) and Magnetically Arrested Disk (MAD; 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 results from these simulations has been obtained by general relativistic ray tracing and radiative transfer algorithms, such as GRay (Chan et al. 2013) and ipole (Mościbrodzka & Gammie 2018).

In this paper, we analyze the GRMHD simulations that were run using harm (see Gammie et al. 2003; Narayan et al. 2012). The simulations include SANE and MAD configurations of the magnetic fields and different black hole spins with corotating or counterrotating 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 the 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

$\begin{eqnarray}&&R={R}_{\mathrm{high}}\displaystyle \frac{{\beta }^{2}}{1+{\beta }^{2}}+\displaystyle \frac{1}{1+{\beta }^{2}}.\end{eqnarray} \tag{ 8 }$

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; and EHT 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}_{\mathrm{BH}}{n}_{e}^{2}$ (see a derivation and discussion in the Appendix).

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 versus 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⁹ M_⊙). For this value of the black hole 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³, amounting to a total of ∼30,000 images in this set. Ray tracing and radiative transfer calculations for this simulation set are carried out using GRay.

(ii) The second set of models (henceforth Set B) explores a comprehensive set of black hole spins and R_high of the accretion plasma, which are used in EHT 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⁹ M_⊙. Each model in this simulation set is used to generate 600–1000 images at a time cadence of 5 GM/c³, amounting to ∼50,000 images in this set. Ray tracing and radiative transfer calculations for this simulation set are 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 effect of the approximation on the lightcurve of flux density in a simulation is restricted to only a few percent. Depending upon the nature of the spatial correlations in structural variability in the simulation, the computed images may carry either time-correlated or time-uncorrelated 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 spaces 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³. 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 (EHT Collaboration et al. 2019e). Gebhardt et al. (2011) measured a mass of M_BH = 6.6 ±0.4 × 10⁹ M_⊙ using stellar dynamics. Walsh et al. (2013) measured a mass of ${M}_{\mathrm{BH}}={3.5}_{-0.7}^{+0.9}\times {10}^{9}{M}_{\odot }$ 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. (2011; see EHT Collaboration et al. 2019f). Because of this, we allow the black hole mass to vary in the range 5 × 10⁹ M_⊙ ≤M_BH ≤ 9.5 × 10⁹ 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 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 EHT 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⁴ GM/c³. 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³ ≤ t ≤ 10⁴ GM/c³.

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 two-dimensional Fourier transform on each of the images. We 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

$\begin{eqnarray}&&D=\displaystyle \frac{1}{n}\sum _{i=1}^{n}\left\{1-\cos ({\theta }_{i}-\bar{\theta })\right\},\end{eqnarray} \tag{ 9 }$

where $\bar{\theta }$ represents the circular mean, defined as

$\begin{eqnarray}&&\bar{\theta }={\tan }^{-1}\left\{\displaystyle \frac{\sum _{i=1}^{n}\sin ({\theta }_{i})}{\sum _{i=1}^{n}\cos ({\theta }_{i})}\right\}.\end{eqnarray} \tag{ 10 }$

In the limit of small deviations in the time-series from the circular mean $\bar{\theta }$, the directional dispersion can be approximated as

$\begin{eqnarray}&&D\approx \displaystyle \frac{{\sigma }^{2}}{2},\end{eqnarray} \tag{ 11 }$

where σ is inferred as a standard deviation of the time-series θ_i for small fluctuations about $\bar{\theta }$.

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

Zoom In Zoom Out Reset image size

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.

4.1. 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 time span 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 time span of Δt = 6 days using Equation (11) as

$\begin{eqnarray}&&{\sigma }^{2}({t}_{0})=\displaystyle \frac{1}{{\rm{\Delta }}t}{\int }_{{t}_{0}}^{{t}_{0}+{\rm{\Delta }}t}{dt}\ 2\left[1-\cos \{{\phi }_{\mathrm{cp}}(t)-{\bar{\phi }}_{\mathrm{cp}}\}\right].\end{eqnarray} \tag{ 12 }$

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 observations, 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 the closure phase with time during one segment of a simulation chosen from Set B, which 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 are comparable to the length of the observations. Given the low degree of variability observed for some triangles in M87 data (see Section 2), this motivates us to define a metric that quantifies how common such periods are in a given model.

Zoom In Zoom Out Reset image size

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 σ₁, σ₂, σ₃ of this observed variability. In Figure 5, we choose as the most likely values of ${\sigma }_{{c}_{0}}$ 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 ${\mathscr{F}}({\sigma }_{1},{\sigma }_{2},{\sigma }_{3})$ of a model as the fraction of 6 day segments that shows 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.

Zoom In Zoom Out Reset image size

Formally, ${\mathscr{F}}({\sigma }_{1},{\sigma }_{2},{\sigma }_{3})$ can be written as

$\begin{eqnarray}\begin{array}{rcl}{\mathscr{F}}({\sigma }_{1},{\sigma }_{2},{\sigma }_{3}) & = & {\displaystyle \int }_{0}^{\infty }d{\sigma }_{1}^{{\prime} }{\displaystyle \int }_{0}^{\infty }d{\sigma }_{2}^{{\prime} }{\displaystyle \int }_{0}^{\infty }d{\sigma }_{3}^{{\prime} }\ P({\sigma }_{\mathrm{sys}})\\ & & \times \,{P}_{\mathrm{model}}({\sigma }_{1}^{{\prime} },{\sigma }_{2}^{{\prime} },{\sigma }_{3}^{{\prime} }| {\sigma }_{\mathrm{sys}}),\end{array}\end{eqnarray} \tag{ 13 }$

where P(σ_sys) is the prior in the systematic uncertainty σ_sys, and ${P}_{\mathrm{model}}({\sigma }_{1}^{{\prime} },{\sigma }_{2}^{{\prime} },{\sigma }_{3}^{{\prime} }| {\sigma }_{\mathrm{sys}})$ is the three-dimensional distribution of circular standard deviations ${\sigma }_{1}^{{\prime} },{\sigma }_{2}^{{\prime} },{\sigma }_{3}^{{\prime} }$ computed in the three low-variability triangles using Equation (12). The observed variability, as quantified in Section 2, is a combination of the true intrinsic variability of the source and the systematic errors in the 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

$\begin{eqnarray}\begin{array}{rcl}{\mathscr{F}}({\sigma }_{1},{\sigma }_{2},{\sigma }_{3}) & = & {\int }_{0}^{{\sigma }_{1}}d{\sigma }_{1}^{{\prime} }{\int }_{0}^{{\sigma }_{2}}d{\sigma }_{2}^{{\prime} }\\ & & \times \,{\int }_{0}^{{\sigma }_{3}}d{\sigma }_{3}^{{\prime} }\ {P}_{\mathrm{model}}({\sigma }_{1}^{{\prime} },{\sigma }_{2}^{{\prime} },{\sigma }_{3}^{{\prime} }),\end{array}\end{eqnarray} \tag{ 14 }$

i.e., simply assuming that the combined intrinsic degree of variability and the systematic uncertainties cannot exceed the inferred values.

Because our inferred degree of variability from the observations themselves have associated formal errors as described by the marginalized posteriors P(σ₁), P(σ₂), and P(σ₃) as shown in Figure 3; we compute the compliance fraction by integrating over those posteriors using

$\begin{eqnarray}&&{\mathscr{C}}(\mathrm{model})=\int d{\sigma }_{1}d{\sigma }_{2}d{\sigma }_{3}P({\sigma }_{1})P({\sigma }_{2})P({\sigma }_{3}){\mathscr{F}}({\sigma }_{1},{\sigma }_{2},{\sigma }_{3}).\end{eqnarray} \tag{ 15 }$

The compliance fraction denotes the expectation value of ${\mathscr{F}}({\sigma }_{1},{\sigma }_{2},{\sigma }_{3})$ given the probability distributions of σ₁, σ₂, and σ₃.

We compute this compliance fraction for black hole masses in the range 5 × 10⁹ M_⊙ ≤ M_BH ≤ 9.5 × 10⁹ M_⊙ and use the maximum number as the compliance fraction for a given GRMHD model.

In this section, we use the two sets of simulations we discuss in Section 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 the closure phase across the six 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 a smaller difference in compliance fractions between the MAD and SANE simulations. We observe that, for small values of R_high, the dependence of the 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 EHT Collaboration et al. (2019d, 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 Figure 7 (Chan et al. 2015).

Zoom In Zoom Out Reset image size

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 across the six days observed in M87 argues in favor of images characterized by thin (∼20%) rings of emission.

Zoom In Zoom Out Reset image size

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 corotating, 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 (EHT 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⁴⁰ erg s⁻¹, 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⁴² to 10⁴⁵ erg s⁻¹. Simulations with fast black hole spins, which are consistent with these constraints, are also characterized by larger compliance fractions, giving preference to these models.

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 observation 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 (EHT 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 the 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 the black hole spin, with the 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 the day-to-day variability in the closure-phase data and compare it to the models is a useful tool in exploring the origin of variability in horizon-scale images of black holes and in discriminating between models.

This work was supported, in part, by the NSF PIRE award 1743747 and NASA ATP award 80NSSC20K0521. L.M. acknowledges support from an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award No. AST-1903847. M.W. acknowledges the support of the Black Hole Initiative at Harvard University, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation to Harvard University. All ray tracing calculations for Set A were performed with the El Gato GPU cluster at the University of Arizona that is funded by NSF award 1228509. All analyses for Set B were performed on CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442.

The EHT Collaboration thanks the following organizations and programs: the Academy of Finland (projects 274477, 284495, 312496, and 315721); the Agencia Nacional de Investigación y Desarrollo, Chile via NCN19₀₅₈ (TITANs) and Fondecyt 3190878; the Alexander von Humboldt Stiftung; an Alfred P. Sloan Research Fellowship; Allegro, the European ALMA Regional Centre node in the Netherlands, the NL astronomy research network NOVA, and the astronomy institutes of the University of Amsterdam, Leiden University, and Radboud University; the black hole Initiative at Harvard University, through a grant (60477) from the John Templeton Foundation; the China Scholarship Council; Consejo Nacional de Ciencia y Tecnología (Mexico, projects U0004-246083, U0004-259839, F0003-272050, M0037-279006, F0003-281692, 104497, 275201, and 263356); the Delaney Family via the Delaney Family John A. Wheeler Chair at Perimeter Institute; Dirección General de Asuntos del Personal Académico-Universidad Nacional Autónoma de México (projects IN112417 and IN112820); the European Research Council Synergy Grant "BlackHoleCam: Imaging the Event Horizon of Black Holes" (grant 610058); the Generalitat Valenciana postdoctoral grant APOSTD/2018/177 and GenT Program (project CIDEGENT/2018/021); MICINN Research Project PID2019-108995GB-C22; the Gordon and Betty Moore Foundation (grant GBMF-3561); the Istituto Nazionale di Fisica Nucleare sezione di Napoli, iniziative specifiche TEONGRAV; the International Max Planck Research School for Astronomy and Astrophysics at the Universities of Bonn and Cologne; Joint Princeton/Flatiron and Joint Columbia/Flatiron Postdoctoral Fellowships; research at the Flatiron Institute is supported by the Simons Foundation; the Japanese Government (Monbukagakusho: MEXT) Scholarship; the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for JSPS Research Fellowship (JP17J08829); the Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS; grants QYZDJ-SSW-SLH057, QYZDJSSW-SYS008, and ZDBS-LY-SLH011); the Leverhulme Trust Early Career Research Fellowship; the Max-Planck-Gesellschaft (MPG); the Max-Planck Partner Group of the MPG and the CAS; the MEXT/JSPS KAKENHI (grants 18KK0090, JP18K13594, JP18K03656, JP18H03721, 18K03709, 18H01245, and 25120007); the Malaysian Fundamental Research Grant Scheme (FRGS). FRGS/1/2019/STG02/UM/02/6; the MIT International Science and Technology Initiatives Funds; the Ministry of Science and Technology (MOST) of Taiwan (105-2112-M-001-025-MY3, 106-2112-M-001-011, 106-2119- M-001-027, 107-2119-M-001-017, 107-2119-M-001-020, 107-2119-M-110-005, 108-2112-M-001-048, and 109-2124-M-001-005); the National Aeronautics and Space Administration (NASA; Fermi Guest Investigator grant 80NSSC20K1567, NASA Astrophysics Theory Program grant 80NSSC20K0527, and NASA NuSTAR award 80NSSC20K0645); the National Institute of Natural Sciences of Japan; the National Key Research and Development Program of China (grants 2016YFA0400704 and 2016YFA0400702); the National Science Foundation (NSF; grants AST-0096454, AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-1126433, AST-1140030, DGE-1144085, AST-1207704, AST-1207730, AST-1207752, MRI-1228509, OPP-1248097, AST-1310896, AST-1555365, AST-1615796, AST-1715061, AST-1716327, AST-1903847, and AST-2034306); the Natural Science Foundation of China (grants 11573051, 11633006, 11650110427, 10625314, 11721303, 11725312, 11933007, 11991052, and 11991053); a fellowship of China Postdoctoral Science Foundation (2020M671266); the Natural Sciences and Engineering Research Council of Canada (NSERC; including a Discovery Grant and the NSERC Alexander Graham Bell Canada Graduate Scholarships-Doctoral Program); the National Youth Thousand Talents Program of China; the National Research Foundation of Korea (the Global PhD Fellowship Grant: grants NRF-2015H1A2A1033752 and 2015- R1D1A1A01056807; the Korea Research Fellowship Program: NRF-2015H1D3A1066561 and Basic Research Support Grant 2019R1F1A1059721); the Netherlands Organization for Scientific Research VICI award (grant 639.043.513) and Spinoza Prize SPI 78-409; the New Scientific Frontiers with Precision Radio Interferometry Fellowship awarded by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Technology of South Africa; the Onsala Space Observatory (OSO) national infrastructure, for the provisioning of its facilities/observational support (OSO receives funding through the Swedish Research Council under grant 2017–00648) the Perimeter Institute for Theoretical Physics (research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science, and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation, and Science); the Spanish Ministerio de Economía y Competitividad (grants PGC2018-098915-B-C21, AYA2016-80889-P, and PID2019-108995GB-C21); the State Agency for Research of the Spanish MCIU through the "Center of Excellence Severo Ochoa" award for the Instituto de Astrofísica de Andalucía (SEV-2017–0709); the Toray Science Foundation; the Consejería de Economía, Conocimiento, Empresas y Universidad of the Junta de Andalucía (grant P18-FR-1769), the Consejo Superior de Investigaciones Científicas (grant 2019AEP112); the US Department of Energy (USDOE) through the Los Alamos National Laboratory (operated by Triad National Security, LLC), for the National Nuclear Security Administration of the USDOE (Contract 89233218CNA000001); the European Union's Horizon 2020 research and innovation program under grant agreement No. 730562 RadioNet; ALMA North America Development Fund; the Academia Sinica; Chandra DD7-18089X and TM6-17006X; the GenT Program (Generalitat Valenciana) Project CIDEGENT/2018/021. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562, and CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442. XSEDE Stampede2 resource at TACC was allocated through TG-AST170024 and TG-AST080026N. XSEDE JetStream resource at PTI and TACC was allocated through AST170028. The simulations were performed in part on the SuperMUC cluster at the LRZ in Garching, on the LOEWE cluster in CSC in Frankfurt, and on the HazelHen cluster at the HLRS in Stuttgart. This research was enabled, in part, by support provided by Compute Ontario (http://computeontario.ca), Calcul Quebec (http://www.calculquebec.ca), and Compute Canada (http://www.computecanada.ca). We thank the staff at the participating observatories, correlation centers, and institutions for their enthusiastic support. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.01154.V. ALMA is a partnership of the European Southern Observatory (ESO; Europe, representing its member states), NSF, and the National Institutes of Natural Sciences of Japan, together with National Research Council (Canada), MOST (Taiwan), the Academia Sinica Institute of Astronomy and Astrophysics (ASIAA; Taiwan), and the Korea Astronomy and Space Science Institute (KASI; Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, Associated Universities, Inc. (AUI)/NRAO, and the National Astronomical Observatory of Japan (NAOJ). The NRAO is a facility of the NSF operated under cooperative agreement by AUI. APEX is a collaboration between the Max-Planck-Institut für Radioastronomie (Germany), ESO, and the OSO (Sweden). The SMA is a joint project between the SAO and ASIAA and is funded by the Smithsonian Institution and the Academia Sinica. The JCMT is operated by the East Asian Observatory on behalf of the NAOJ, ASIAA, and KASI, as well as the Ministry of Finance of China, CAS, and the National Key R&D Program (No. 2017YFA0402700) of China. Additional funding support for the JCMT is provided by the Science and Technologies Facility Council (UK) and participating universities in the UK and Canada. The LMT is a project operated by the Instituto Nacional de Astrófisica, Óptica, y Electrónica (Mexico) and the University of Massachusetts at Amherst (USA). The IRAM 30 m telescope in Pico Veleta, Spain is operated by IRAM and supported by CNRS (Centre National de la Recherche Scientifique, France), MPG (Max-Planck Gesellschaft, Germany), and IGN (Instituto Geográfico Nacional, Spain). The SMT is operated by the Arizona Radio Observatory, a part of the Steward Observatory of the University of Arizona, with financial support of operations from the State of Arizona and financial support for instrumentation development from the NSF. The SPT is supported by the NSF through grant PLR-1248097. Partial support is also provided by the NSF Physics Frontier Center grant PHY-1125897 to the Kavli Institute of Cosmological Physics at the University of Chicago, the Kavli Foundation, and the Gordon and Betty Moore Foundation grant GBMF 947. The SPT hydrogen maser was provided on loan from the GLT, courtesy of ASIAA. The EHTC has received generous donations of FPGA chips from Xilinx Inc., under the Xilinx University Program. The EHTC has benefited from technology shared under an open-source license by the Collaboration for Astronomy Signal Processing and Electronics Research. The EHT project is grateful to T4Science and Microsemi for their assistance with hydrogen masers. This research has made use of NASA's Astrophysics Data System. We gratefully acknowledge the support provided by the extended staff of the ALMA, both from the inception of the ALMA Phasing Project through the observational campaigns of 2017 and 2018. We would like to thank A. Deller and W. Brisken for the EHT-specific support with the use of DiFX. We acknowledge the significance that Maunakea, where the SMA and JCMT-EHT stations are located, has for the indigenous Hawaiian people.

In this appendix, we present analytical arguments to show that the mm-wavelength images calculated from GRMHD simulations of accretion flows with parameters relevant to the M87 black hole are approximately invariant to the product ${n}_{e}^{2}{M}_{\mathrm{BH}}$, where n_e is the electron-density scale and M_BH is the mass of the black hole.

The radiative transfer equation at a frequency ν is given by

$\begin{eqnarray}&&\displaystyle \frac{{{dI}}_{\nu }}{{ds}}={\eta }_{\nu }-{\chi }_{\nu }{I}_{\nu },\end{eqnarray} \tag{ A1 }$

where I_ν denotes the specific intensity at a frequency ν, s denotes the distance traveled along the path of a photon, η_ν denotes the emissivity of the medium, and χ_ν denotes the opacity of the medium. Equation (A1) can be recast as

$\begin{eqnarray}&&\displaystyle \frac{{{dI}}_{\nu }}{{ds}}={\chi }_{\nu }\left(\displaystyle \frac{{\eta }_{\nu }}{{\chi }_{\nu }}-{I}_{\nu }\right).\end{eqnarray} \tag{ A2 }$

The fraction η_ν /χ_ν is the source function S_ν and, for thermal emission, is equal to the blackbody function, which depends solely on temperature T. We denote this as ${\eta }_{\nu }/{\chi }_{\nu }\equiv {{\mathscr{B}}}_{\nu }(T)$.

Throughout this work, we use the analytic approximation devised by Leung et al. (2011) for synchrotron emissivity from a thermal distribution of relativistic electrons

$\begin{eqnarray}&&{\eta }_{\nu }={n}_{e}\displaystyle \frac{\sqrt{2}\pi {e}^{2}{\nu }_{{\rm{s}}}}{3{K}_{2}(1/{{\rm{\Theta }}}_{e})c}{\left({X}^{1/2}+{2}^{11/12}{X}^{1/6}\right)}^{2}\exp (-{X}^{1/3}),\end{eqnarray} \tag{ A3 }$

where X ≡ ν/ν_s, ${\nu }_{{\rm{s}}}\equiv (2/9){\nu }_{{\rm{c}}}{{\rm{\Theta }}}_{e}^{2}\sin \theta$, ν_c ≡ eB/(2π m_e c), Θ_e ≡k_b T_e/(m_e c²) is the dimensionless electron temperature, and K₂ denotes the modified Bessel function of the second kind. In the these expressions, θ denotes the angle between the magnetic field and the emitted photon, B denotes the magnetic field strength, and T_e denotes the electron temperature. We use the constants m_e , e, k_b, and c to denote the electron mass, electron charge, the Boltzmann's constant, and the speed of light respectively. It is explicit in this expression that the synchrotron emissivity and opacity, for a given electron temperature and magnetic field, scale proportionally to the electron number density n_e. We show here that, for the parameters of the black hole at the center of the M87 galaxy, the synchrotron opacity and emissivity at a wavelength of λ = 1.3 mm also scale proportionally to a power of the magnetic field, i.e., ∝ B^α , with α ∼ 2.

We estimate the properties of the plasma in the inner accretion flow of M87 using the following analytic model. The electron density at a given radius r is given by the continuity equation, assuming a spherical accretion rate $\dot{M}$ via the relation

$\begin{eqnarray}&&\dot{M}=4\pi {r}^{2}\left(\displaystyle \frac{h}{r}\right){m}_{p}{n}_{e}{u}^{r},\end{eqnarray} \tag{ A4 }$

where h/r is the scale height of the accretion flow, m_p is the proton mass (assuming fully ionized hydrogen), and u^r is the radial component of the accretion flow. We take the latter to be a fraction ξ of the freefall velocity, i.e.,

$\begin{eqnarray}&&{u}^{r}=\xi {\left(\displaystyle \frac{{{GM}}_{\mathrm{BH}}}{r}\right)}^{1/2}.\end{eqnarray} \tag{ A5 }$

We also express the accretion rate as $\dot{M}\equiv \dot{m}{\dot{M}}_{{\rm{E}}}$, in terms of the Eddington accretion rate defined by

$\begin{eqnarray}&&{\dot{M}}_{{\rm{E}}}\equiv \displaystyle \frac{{L}_{{\rm{E}}}}{\epsilon {c}^{2}}=\displaystyle \frac{4\pi {{GM}}_{\mathrm{BH}}{m}_{{\rm{p}}}}{\epsilon c{\sigma }_{{\rm{T}}}},\end{eqnarray} \tag{ A6 }$

where L_E is the Eddington critical luminosity, is the radiative efficiency of the flow, σ_T is the Thomson cross section, and $\dot{m}$ is the mass accretion rate in units of the Eddington accretion rate. Under these conditions, the electron density is

$\begin{eqnarray}&&\begin{array}{l}{n}_{e}=\left(\displaystyle \frac{{c}^{2}}{{{GM}}_{\mathrm{BH}}{\sigma }_{{\rm{T}}}}\right)\dot{m}{\epsilon }^{-1}{\xi }^{-1}{\left(\displaystyle \frac{h}{r}\right)}^{-1}{\left(\displaystyle \frac{{{rc}}^{2}}{{{GM}}_{\mathrm{BH}}}\right)}^{-3/2}\\ =\,7\times {10}^{5}{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{6.5\times {10}^{9}{M}_{\odot }}\right)}^{-1}\left(\displaystyle \frac{\dot{m}}{2\times {10}^{-5}}\right)\\ \ {\left(\displaystyle \frac{{{rc}}^{2}}{5{{GM}}_{\mathrm{BH}}}\right)}^{-3/2}\,{\mathrm{cm}}^{-3},\end{array}\end{eqnarray} \tag{ A7 }$

where in the last expression we have used typical values for the M87 black hole and set ξ = = 0.1 and h/r = 0.4.

Because of the radiatively inefficient character of the accretion flow around the M87 black hole, the ion temperature T_i is a fraction of the virial temperature

$\begin{eqnarray}&&{T}_{{\rm{v}}}=\displaystyle \frac{{{GM}}_{\mathrm{BH}}{m}_{p}}{3{k}_{{\rm{B}}}r}.\end{eqnarray} \tag{ A8 }$

It is also customary to write the electron temperature as a fraction 1/R of the ion temperature, i.e.,

$\begin{eqnarray}&&{T}_{{\rm{e}}}=\displaystyle \frac{1}{R}\left(\displaystyle \frac{{{GM}}_{\mathrm{BH}}{m}_{{\rm{p}}}}{3{k}_{{\rm{B}}}r}\right)\left(\displaystyle \frac{{T}_{i}}{{T}_{{\rm{v}}}}\right).\end{eqnarray} \tag{ A9 }$

Finally, we write the magnetic field in terms of the plasma β parameter as

$\begin{eqnarray}&&{n}_{e}{k}_{{\rm{B}}}({T}_{i}+{T}_{e})=\beta \displaystyle \frac{{B}^{2}}{8\pi }\end{eqnarray} \tag{ A10 }$

such that

$\begin{eqnarray}&&\begin{array}{l}B=2{c}^{2}{\left(\displaystyle \frac{2\pi {m}_{{\rm{p}}}}{3{{GM}}_{\mathrm{BH}}{\sigma }_{{\rm{T}}}}\right)}^{1/2}{\left[\displaystyle \frac{\dot{m}(1+R)}{\epsilon \xi \beta R}\right]}^{1/2}\\ \ \times \,{\left(\displaystyle \frac{{T}_{i}}{{T}_{{\rm{v}}}}\right)}^{1/2}{\left(\displaystyle \frac{h}{r}\right)}^{-1/2}{\left(\displaystyle \frac{{{rc}}^{2}}{{{GM}}_{\mathrm{BH}}}\right)}^{-5/4}\\ \ =\,9{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{6.5\times {10}^{9}{M}_{\odot }}\right)}^{-1/2}{\left(\displaystyle \frac{\dot{m}}{2\times {10}^{-5}}\right)}^{1/2}{\left(\displaystyle \frac{{{rc}}^{2}}{5{{GM}}_{\mathrm{BH}}}\right)}^{-5/4}\,{\rm{G}},\end{array}\end{eqnarray} \tag{ A11 }$

and we have set T_i /T_v = 1/3 and β = 10.

Figure 9 shows the ratio η_ν /n_e evaluated at the observed frequency of the EHT (ν = 230 GHz, λ = 1.33 mm) and at a radial distance r = 5GM_BH/c², as a function of the strength of the magnetic field B, for the values of the various parameters used in the previous equations, i.e., for the conditions of the inner accretion flow around the black hole at the center of M87. This figure demonstrates that, over a broad range of magnetic field strengths (spanning more than an order of magnitude), the synchrotron emissivity scales approximately as η_ν ∼ n_e B^α , with α ≃ 2. This allows us to write the emissivity as

$\begin{eqnarray}&&{\eta }_{\nu }\simeq {n}_{{\rm{e}}}{B}^{\alpha }{{ \mathcal B }}_{\nu }(T){f}_{\nu }(T,{\boldsymbol{r}}),\end{eqnarray} \tag{ A12 }$

and the opacity as

$\begin{eqnarray}&&{\chi }_{\nu }\simeq {n}_{{\rm{e}}}{B}^{\alpha }{f}_{\nu }(T,{\boldsymbol{r}}),\end{eqnarray} \tag{ A13 }$

where f_ν (T, r ) captures the scaling of the opacity with temperature T and the location in the accretion flow (see similar arguments in Chan et al. 2015).

Zoom In Zoom Out Reset image size

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 ${ \mathcal M }$, as long as the magnetic field is also rescaled by a factor of ${{ \mathcal M }}^{1/2}$ and the internal energy by a factor of ${ \mathcal M }$. In other words, $B\sim {n}_{{\rm{e}}}^{1/2}$ and the Alfvén speed $B/{({m}_{{\rm{p}}}{n}_{e})}^{1/2}$ are 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 simulation outputs are invariant to rescaling, as long as the product ${n}_{e}{B}^{\alpha }\sim {n}_{e}^{1+\alpha /2}$ 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.,

$\begin{eqnarray}&&{ds}\equiv {ds}^{\prime} \left(\displaystyle \frac{{{GM}}_{\mathrm{BH}}}{{c}^{2}}\right).\end{eqnarray} \tag{ A14 }$

Rewriting Equation (A2) using the scaling (A12)–(A14), we find

$\begin{eqnarray}&&\displaystyle \frac{{dI}}{{ds}^{\prime} }={n}_{{\rm{e}}}^{1+\alpha /2}f(T,{\boldsymbol{r}})\displaystyle \frac{{{GM}}_{\mathrm{BH}}}{{c}^{2}}\left[{\mathscr{B}}(T)-I\right]\end{eqnarray} \tag{ A15 }$

$\begin{eqnarray}&&=\left({n}_{{\rm{e}}}^{1+\alpha /2}{M}_{\mathrm{BH}}\right)\displaystyle \frac{G}{{c}^{2}}f(T,{\boldsymbol{r}})\left[{\mathscr{B}}(T)-I\right].\end{eqnarray} \tag{ A16 }$

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}_{{\rm{e}}}^{1+\alpha /2}{M}_{\mathrm{BH}}$ and with α ≃ 2.