Future Prospects for Constraining Black Hole Spacetime: Horizon-scale Variability of Astrophysical Jets

The Event Horizon Telescope (EHT) Collaboration has recently published the first horizon-scale images of the supermassive black holes M87* and Sgr A* and provided some first information on the physical conditions in their vicinity. The comparison between the observations and the three-dimensional general relativistic magnetohydrodynamic (GRMHD) simulations has enabled the EHT to set initial constraints on the properties of these black hole spacetimes. However, accurately distinguishing the properties of the accretion flow from those of the spacetime, most notably, the black hole mass and spin, remains challenging because of the degeneracies the emitted radiation suffers when varying the properties of the plasma and those of the spacetime. The next-generation EHT (ngEHT) observations are expected to remove some of these degeneracies by exploring the complex interplay between the disk–jet dynamics, which represents one of the most promising tools for extracting information on the black hole spin. By using GRMHD simulations of magnetically arrested disks and general relativistic radiative transfer (GRRT) calculations of the emitted radiation, we have studied the properties of the jet and the accretion disk dynamics on spatial scales that are comparable with the horizon. In this way, we are able to highlight that the radial and azimuthal dynamics of the jet are well correlated with the black hole spin. Based on the resolution and image reconstruction capabilities of the ngEHT observations of M87*, we can assess the detectability and associated uncertainty of this correlation. Overall, our results serve to assess the prospects for constraining the black hole spin with future EHT observations.


INTRODUCTION
The observational exploration of the physical conditions in the vicinity of a black hole is one of the most promising avenues for an accurate mapping of the spacetime and understanding of the origin of the relativistic jets.Through groundbreaking efforts, the Event Horizon Telescope (EHT) Collaboration has recently provided images of the horizon-scale emission of the accretion flow in the vicinity of supermassive black holes in the nearby radio galaxy M87 (M87 * ) and in the Galactic center (Sgr A * ).More specifically, in 2017, the EHT array consisted of eight stations at six geographical sites that successfully observed M87 * (Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f) and Sgr A * (Event Horizon Telescope Collaboration et al. 2022a,b,c,d,e,f) with the unprecedented horizon-scale resolution of ∼ 20 µas.The imaging surveys of M87 * and Sgr A * revealed bright emission rings that are consistent with the theoretical predictions of general relativity.
Comparing the EHT observations with the theoretical simulations enables us to constrain the black-hole spacetime within the mass and spin with general relativity as well as alternative theories.The measured shadow diameters (d ∼ 42 µas for M87 * and d ∼ 52 µas for Sgr A * ) are consistent with the prediction of general relativity (e.g., Event Horizon Telescope Collaboration et al. 2019e, 2022e;Younsi et al. 2023) and can be also used to discriminate from general relativity's possible deviations of the black-hole spacetime (e.g., Kocherlakota & Rezzolla 2020;Event Horizon Telescope Collaboration et al. 2019f, 2022f;Cruz-Osorio et al. 2021;Vagnozzi et al. 2022).The comparison of the observations with a large library of three-dimensional general-relativistic magnetohydrodynamic (GRMHD) simulations provides information on the plausible ranges of model parameters for the spacetime, the magnetic field, and the accretion flow properties (e.g., Event Horizon Telescope Collaboration et al. 2019e, 2022e).
Furthermore, the horizon-scale detection of the polarized profile around M87 * has provided additional constraints on the magnetic fields (Event Horizon Telescope Collaboration et al. 2021a,b).In particular, the polarization properties of M87 * can be explained with magnetic fields around the emission region that are predominantly poloidal, and this has allowed us to obtain the first measurements of the magnetic field strength, electron density, and temperature with horizon-scale resolution.Moreover, by comparing the imaging results with the library of GRMHD polarized images, it was possible to confirm their consistency with magnetically arrested accretion disks (MADs).
Despite this considerable progress, with the current angular resolution of the EHT observations, it remains challenging to constrain the detailed structure of the black-hole spacetime solely from the horizon-scale images.Despite the sig-nificant progress achieved, the current angular resolution of the EHT observations poses challenges in accurately determining the detailed structure of black-hole spacetime based solely on horizon-scale images.For instance, according to general relativity, the diameter of the black-hole shadow is predicted to depend weakly on the observer's inclination angle and the black-hole spin, approximately 5.0±0.2rg , where r g represents the gravitational radius (e.g., see Takahashi 2004;Event Horizon Telescope Collaboration et al. 2022f).Although these parameters are crucial in characterizing the spacetime of a black hole, achieving a more precise constraint on the black-hole parameters would require additional criteria and/or improvements in observations.The essential role of non-thermal emission has been demonstrated with observations of M87 * at the regime between low (43 and 86 GHz) and high (near-IR) frequencies (Pandya et al. 2016;Davelaar et al. 2019;Cruz-Osorio et al. 2022;Fromm et al. 2022).Overall, GRMHD simulations have clearly demonstrated the intricate dependence of the total intensity and polarization map on the properties of both the plasma and spacetime in a complex manner.Differentiating between them is not a straightforward endeavor and necessitates the application of innovative methodologies.
Future developments of the EHT capabilities as those explored with the next generation EHT (ngEHT) project, represent one of the most promising to address the observational and theoretical challenges discussed above.The aim of the development is to achieve improved angular resolution across a wide range of frequencies, including 86, 230, and 345 GHz.The possibility of improved spacetime measurements utilizing ngEHT observations is currently being intensively explored (e.g., Chael et al. 2021;Ricarte et al. 2023), while some research has investigated the non-thermal emission effects and constrained the observational results.In particular, Fromm et al. (2022) have parameterized the nonthermal effect and investigated the characteristics of model images and broadband spectra spanning from the radio to the near-infrared range for M87 * .The obtained best-bet models for each spin showed a strong agreement with the spectral features observed at multiple wavelengths (Davelaar et al. 2019).
In addition, the horizon-scale imaging of the M87 * jet represents a promising tool for extracting information on the spacetime properties and measuring the non-thermal effects simultaneously.The observation of the M87 * jet has a long history and has been conducted across a broad frequency range, from radio to γ-ray wavelengths (e.g., Curtis 1918;Reid et al. 1982;Hada et al. 2013;Kim et al. 2018;Snios et al. 2019;MAGIC Collaboration et al. 2020).The time-averaged jet morphology was further compared with 86 GHz GMVA observations of M87 * to constrain the plausible parameter spaces of GRMHD and general-relativistic radiative transfer (GRRT) simulations (Cruz-Osorio et al. 2022;Fromm et al. 2022).Especially, recent 86 GHz imaging results of M87 (Lu et al. 2023) have revealed an edgebrightened jet that is connected to the accretion flow of the black hole.This observational milestone provides strong motivation for theoretical and observational studies aimed at investigating the relationship between the jet and photon ring.The recent observational improvement and forthcoming ngEHT observations hold great potential to detect the horizon-scale structure of the jet and to capture the disk dynamics of the M87 * through long-time monitoring (Roelofs et al. 2023).
In this paper, we investigate the dynamics of jets and accretion disks that can be probed through the ngEHT observations.The previous studies in this area (e.g., Raymond et al. 2021;Roelofs et al. 2023) have provided evidence for the detectability of the jet base and photon ring.We report on three-dimensional GRMHD simulations of MAD disks with the numerical code BHAC to obtain an accurate modeling of the accretion flow in the vicinity of the black hole (Porth et al. 2017;Olivares et al. 2019).In a comparative study with other GRMHD codes utilized in the EHT Collaboration, BHAC demonstrated maturity, capability, and consistency in simulation results (Porth et al. 2019).Additionally, a comprehensive convergence study examining the Riemann problem, shocks, wave propagation, and accretion is available in Porth et al. (2017) and Olivares et al. (2019).The synthetic images of the emission are calculated with the GRRT scheme implemented in the numerical code BHOSS, introducing a hybrid thermal-non-thermal particle distribution function optimized by the observed continuum spectra of M87 * (Fromm et al. 2022;Cruz-Osorio et al. 2022).The Radiative Transfer Scheme was rigorously tested to ensure its consistency and accuracy within the context of the general relativistic radiative transfer codes used by the EHT Collaboration (Gold et al. 2020).We primarily focus on the dynamic properties of the jet and accretion disk and demonstrate that the azimuthal variation and radial jet propagation provide information on the black-hole spacetime.Furthermore, by comparing the simulations and reconstructed images, we measure the uncertainty and detectability of the black-hole characteristics under the planned ngEHT observations.The plan of this paper is as follows: In Section 2, we introduce the setup of the GRMHD and GRRT simulations and the strategy of the synthetic ngEHT observations.In Section 3, we demonstrate the fundamental properties of the utilized GRMHD and GRRT models.In Section 4, we investigate the properties of the jet-disk dynamics and demonstrate the detectability with the expected ngEHT observations.Section 5 is instead devoted to investigating the physical origin of the detected jet dynamics discussed in Section 4. Finally, in Section 6, we present our conclusions and future prospects.

General-Relativistic magnetohydrodynamic simulations and radiative transfer
Our investigation of the dynamics of jet launching and accretion flow around a black hole has employed threedimensional GRMHD simulations performed with the numerical code BHAC (Porth et al. 2017;Olivares et al. 2019), which solves the GRMHD equations on a Kerr background expressed in the spherical Kerr-Schild coordinates (t, r, θ, ϕ).The grid spacing is logarithmic in the radial and linear in the polar and azimuthal directions.Hereafter, we use geometrized units, in which the gravitational constant, G, and the speed of light, c, are set to be unity and define the gravitational radius as r g := GM/c 2 , where M is the black-hole mass and the gravitational time t g := r g /c.
The conditions of GRMHD simulations (the convention, resolution of the grid, and initial conditions employed in our analysis) are the same as those in previous simulations, which were used to compare with M87 * observations (Fromm et al. 2022;Cruz-Osorio et al. 2022).We take the time window with a quasi-stationary mass accretion rate between 13000 ≤ t/t g ≤ 15000, with a time resolution of t = 10 t g ∼ 3.6 days for M87 * .The scaling factors of the accretion rate (summarized in Table 1) and model parameters are selected to reproduce a consistent total flux at 230 GHz (∼ 1 Jy) and a continuum spectrum density within the frequency range of 10 10 Hz ≤ ν ≤ 10 16 Hz.
Utilizing the GRMHD results, we calculated synthetic images using the general-relativistic ray-tracing code employed in the BHOSS (Younsi et al. 2012(Younsi et al. , 2020)), which solves the co-variant radiative-transfer equation.The calculation of the synthetic images is performed with the fast light approximation, in which the fluid is assumed not to change during the photon propagation so that an image corresponds to a single time slice.Synchrotron radiation is estimated to be the main (only) source of emission at 230 GHz and the electron energy distribution is characterized by the socalled "kappa" model (Xiao 2006) with a non-thermal energy contribution from the jet (Davelaar et al. 2019).We focus specifically on the best-bet models of MAD with five different spins (a * = −0.9375,−0.5, 0.0, 0.5, and 0.9375, where a * = J/M 2 represents the normalized spin and J denotes the angular momentum) at an inclination angle of i = 160 • .These models exhibit a continuum spectrum consistent with the observations of M87 * (Cruz- Osorio et al. 2022;Event Horizon Telescope Collaboration et al. 2019e).The blackhole mass is set to be 6.5 × 10 9 M ⊙ , and the distance is set to be 16.8 Mpc.Additionally, we set the azimuthal angle of the spin axis as ϕ obs (= arctan(y obs /x obs )) = 18 • (Fromm et al. 2022).As already anticipated, the primary aim of this paper is to explore the dynamics of the jet and of the disk, which contain information about the properties of spacetime, and to assess their detectability with the anticipated capabilities of the ngEHT.In order to comprehensively analyze the detectable properties, it is essential to assess the impact of observational uncertainties with the expected imaging results.The observed visibility V (t, u, v) is defined as the Fourier transform of the source intensity I(t, x obs , y obs ), V (t, u, v) := I(t, x obs , y obs )e 2πi(ux obs +vy obs ) dx obs dy obs , where (x obs , y obs ) are the Cartesian coordinates on the image plane, while (u, v) are the components of the baseline vector (commonly referred to as spatial frequencies) between two antennas, which is normalized by the observing wavelength and projected on the image plane orthogonal to the source direction.The expected telescope properties were examined in detail by Raymond et al. (2021) and Roelofs et al. (2023), and the baseline (u, v) coverage is shown in Fig. 1 for the observational frequency of 230 GHz.The measured visibility on a baseline is the summation of the true visibility and systematic noise, which is assessed with a Gaussian distribution in the complex plane with the standard deviation

Synthetic ngEHT observations and Imaging strategy
where η ≃ 0.88 represents the quantization efficiency (Thompson et al. 2017), ∆ν ≃ 8 GHz refers to the bandwidth of a single spectral window for a specific band and polarization (Roelofs et al. 2023), and ∆t = 10 s denotes the integration time (Event Horizon Telescope Collaboration et al. 2019d).We define the system equivalent flux density (SEFD), as encompassing all thermal noises originating from the telescope receiver chains, Earth's atmosphere, and astronomical background (see Roelofs et al. 2023 for a summary), where the subscripts (i, j) refer to the telescopes involved in the respective baseline.
Assuming that the time duration of a single observation is 24 hours, during which a snapshot image does not change, we generated a data set corresponding to a time window of (2000t g ∼ 720 days).Utilizing the open-source software eht-imaging (Chael et al. 2018a), we generate synthetic data based on the expected ngEHT (u, v) coverage and source models, incorporate systematic noises into the visibilities, and perform calibration on the visibility phases.
We reconstructed the images using the open-source software SMILI, which employs the regularized maximum likelihood (RML) method based on the synthetic data (Akiyama et al. 2017a,b).For the reconstruction, we used visibilities averaged over 6 minutes, setting the imaging field of view to 1200 µas and a pixel size of 4 µas.The reconstruction process was initiated using an image characterized by a circular Gaussian profile with a FWHM of 50 µas.Furthermore, we utilized imaging regularizations of Total Variation (TV), with a parameter value of 10 6 , and relative entropy, with a value of 1 (Event Horizon Telescope Collaboration et al. 2019d).
We present an example of GRRT (top panels) and reconstructed images (middle) in Fig. 2. Each image in the topright and middle-right panels has been convolved with a circular Gaussian function corresponding to the typical resolution of the EHT (20 µas).Note that the reconstructed image accurately captures the extended-jet features, while also exhibiting a certain level of noise structure, which can be attributed to systematic errors and the constrained (u, v) coverage.The intensity ratio between the GRRT and reconstructed images is shown in the bottom panels.The standard deviation of the intensity ratio for the unconvolved images (∼ 2.8) is mitigated through convolution at the nominal resolution (∼ 0.1).In Section 4, we use the convolved images to investigate the common dynamic properties consistent between GRRT and reconstructed images.

Time-averaged image properties
Identifying the spatial regions on the image plane corresponding to the jet and the disk is essential for extracting physical information from both the GRRT and reconstructed Table 2. Time average and standard deviation of jet+disk and disk sizes using the ngEHT sensitivity as defined in Fig. 3, where a * represents the dimensionless black-hole spin, and r obs,,jet+disk and r obs,,disk denote the radial sizes of the jet+disk and disk regions, respectively.images.We calculated the radiation emitted from each region, the boundary of which is defined by the Bernoulli parameter Be := −hu t = 1.02,where h is the specific enthalpy and u t is the time component of the covariant fourvelocity (Fromm et al. 2022).Along each ray trajectory, we set the emissivity to zero outside the corresponding region while maintaining the absorptivity to the relevant value.
To demonstrate the essence of our approach, we show in Fig. 3 the time-averaged properties of the jet and the disk structures as averaged over the entire simulation (i.e., from t/t g = 0 to 2000).The left panel shows the average of the synthetic images for a Kerr black hole with a * = 0.9375 with the green and cyan curves identifying the detectable regions of the disk and jet+disk using the typical ngEHT sensitivity (= 2 × 10 7 K, Doeleman et al. 2019;Raymond et al. 2021;Chael et al. 2021).The disk radiation concentrates within the radial region of the image plane at approximately r obs = x 2 obs + y 2 obs ∼ 100 µas ≃ 26r g , while the jet radiation extends to the radial region 100 µas ≤ r obs < 300 µas.
Each curve in the middle panel of Fig. 3 represents the radial profile of the azimuthally maximum brightness temperature, where the black, blue, and red curves correspond to a * = 0, |a * | = 0.5, and |a * | = 0.9375, respectively.For a * > 0 (solid curve), the radial extent of the jet increases with |a * |: the jet extends towards r obs ∼ 500 µas (300 µas) with a dynamic range of 10 −4 in the spinning (non-spinning) black hole cases (solid curves).In the case of a * < 0 (dashed curves), the radial extent of the jet is similar with a dynamic range of 10 −4 .The existence of a powerful jet in the case of a prograde black hole was also reported by other studies (e.g., Blandford & Znajek 1977;Tchekhovskoy & McKinney 2012;Fromm et al. 2022;Narayan et al. 2022).
The right panel of Fig. 3 shows the radial profile of the time-averaged disk emission.The profile highlights the photon-ring region, with fainter components at the radial center and peak components at 20 µas = 5.2 r g .These features indicate the black-hole shadow and the photon ring, respectively.The disk extends radially from 20 µas to 100 µas, and the brightness temperature within the radial region of r obs ∼ 100 r g decreases.The disk extends radially from ).The green and cyan curves represent the disk and jet+disk regions, respectively, which can be detected with the typical ngEHT sensitivity (= 2 × 10 7 K).In the middle panel, the radial profile of the azimuthal maximum brightness temperature of the jet+disk is shown.The expected sensitivity of the 2017 EHT and ngEHT observations is presented as horizontal gray lines.In the right panel, the same format is followed as in the middle panel, but for the disk emission defined by the Bernoulli parameter (Be < 1.02).The vertical gray lines in the middle and right panels indicate a radius of 100 µas, which is the typical radius of the disk region for each spin case.
20 µas to 100 µas, and the brightness temperature within the radial region of r obs ∼ 100 r g decreases outwards.
The detectable regions of the disk and jet components can be assessed using the sensitivities of the 2017 EHT and the expected ngEHT observations, as represented by the gray horizontal lines in the middle and right panels.The emission from the disk and the jet regions can be detected within the ngEHT sensitivity at radii of r obs < 100 µas and 100 µas ≲ r obs ≲ 300 µas respectively.This would mark a significant improvement compared to the capabilities of the 2017 EHT observations.The uncertainties related to time variation are summarized in Table 2.In light of Fig. 3 and Table 2, we define the radii of the disk and the jet regions on the image plane as r obs < 100 µas and 100 µas ≤ r obs , respectively.

Light-curve properties
We report in Fig. 4 the temporal variability of the emission from all regions, including the jet and the disk.The top panels of Fig. 4 show the light curves for different black-hole spins (a * = −0.9375,−0.5, 0.0, 0.5, and 0.9375), as emitted from all regions (black curve), from the jet (blue) and from the disk (red).The light curves exhibit variability within the entire simulation time window, defined as τ = 2000t g .The mean and standard deviation of the total jet+disk light curve are comparable across each spin case (∼ 0.96 ± 0.16 Jy) and are consistent with the core flux of M87 * observed in previous EHT observations (∼ 1.1−1.2Jy, Doeleman et al. 2012;Event Horizon Telescope Collaboration et al. 2019d).
Table 3. Summary of the quasi-periodic oscillation (QPO).Here, ω represents the peak frequency detected in the power spectral densities (as seen in Fig. 4), Q denotes the quality factor of the QPO, and ∆ω corresponds to the full width at half-maximum (FWHM) of the peak.
Most models demonstrate a dominance of jet radiation over that from the disk.With a ratio of the jet flux to the disk flux that is ⟨F jet ⟩ t /⟨F disk ⟩ t = 2.6, 1.9, 1.1, 0.9, and 1.8 for a * = −0.9375,−0.5, 0.0, 0.5, and 0.9375, respectively, where the brackets ⟨ ⟩ t represent the time average.For both a * ≥ 0.5 and a * < 0 cases, the ratio monotonically increases with the absolute value of the spin, indicating the generation of a powerful jet for rapidly rotating black holes.
Our best-bet GRMHD models have smaller fluxes in the mid-plane region around the black hole (with a field of view of approximately 100 µas and a polar angle range of 57.3 • ≤ θ ≤ 122.7 • ) compared to those of the thermal models based on the EHT GRMHD library (Event Horizon Telescope Collaboration et al. 2019e).To investigate this difference, we calculated the ratio of the total flux in the midplane to that in all regions, utilizing an 80 µas field of view, ).This discrepancy can be attributed to the impact of non-thermal particles, which result in a comparable continuum spectrum spanning from the radio to the near-infrared range and influence the jet morphology at a frequency of 86 GHz in the case of M87 * (see Fromm et al. 2022;Cruz-Osorio et al. 2022).
The variability characteristics of each light curve are also apparent in the power spectral densities (PSDs) P ω , where ω denotes the angular frequency (bottom panels of Fig. 4).The frequency profiles of the PSDs, that are characterized by P ω ∝ ω −2.4±0.5 , are broadly consistent with those produced by red noise (P ω ∝ ω −2 ).This property was already reported in previous research that utilized a large library of GRMHD simulation models without non-thermal particles (Georgiev et al. 2022).
Interestingly, some PSDs display a feature that could be associated with a quasi-periodic oscillation (QPO).This trend is especially pronounced for a * = 0.9375, where the PSD exhibits peak components at time scales of 1/ω ∼ 320t g and 170t g with corresponding quality factors of ω/∆ω ∼ 8 and 2, respectively.Here, ∆ω represents the full width at half maximum (FWHM), determined by fitting the peaks with the Lorentzian function.This pattern is reminiscent of QPO phenomena observed in some stellar-mass and supermassive black holes and could be employed to deduce the properties of the background spacetime (e.g., Abramowicz et al. 2003;Rezzolla et al. 2003;Reynolds & Miller 2009;Kato et al. 2008;McKinney et al. 2012;Smith et al. 2021).We intend to undertake a detailed investigation of the QPOs in future works.

RESULTS: MOVIE MORPHOLOGIES
In what follows we discuss the analysis of the dynamics of the disk and of the jet by employing synthetic images generated from GRRT simulations.Our assessment of the observational detectability and uncertainty has been achieved through the comprehensive imaging analysis of the expected ngEHT observations, which we summarized in Section 2.2.Our goal is to demonstrate that the structures of the jet and of the disk exhibit azimuthal and radial fluctuations that depend on the black-hole spin.Figure 5. Azimuthal profiles of the image brightness radially averaged within the jet region (100 µas ≤ r obs < 300 µas).From left to right, the different panels refer to black holes with a * = −0.9375 to 0.9375.The top and middle panels depict the brightness temperature on the plane of the azimuthal angle and time (ϕ, t), utilizing the GRRT and reconstructed images, respectively.The white points denote the maximum azimuthal brightness at each time.The bottom panels exhibit the distributions of the peak azimuthal angles obtained from the upper panels.The blue and red histograms represent the estimates derived from the GRRT and reconstructed images, respectively.

Variability of the Azimuthal Distribution
In Fig. 5, we show the azimuthal profiles of the jet, where each image is radially averaged within the jet region; i.e., 100 µas ≤ r obs < 300 µas (corresponding to 26 ≤ r obs /r g < 79).The azimuthal profile corresponding to a negative spin exhibits substantial variability along the temporal direction, while profiles with positive spins display relatively stable behavior over time.The white points ϕ obs, peak (t) illustrate the temporal evolution of the azimuthal angle at maximum brightness.In the instance of a * = −0.9375,intermittent azimuthal rotations of ϕ obs, peak (t) occur in the same direction as the accretion flow (clockwise) with a yearly timescale (∼ 1000t g ≃ 1.0 year for M = 6.5×10 9 M ⊙ ).For example, the white points within the range of −180 • ≤ ϕ obs ≤ 60 • and 170 ≤ t/t g ≤ 880 exhibit an angular velocity of dϕ obs /dt ∼ −0.4 • /t g .The substantial variation in the azimuthal angle on a yearly scale, δϕ obs ∼ 240 • , is a distinctive characteristic of the case with a * = −0.9375,separating it from the jets produced by black holes with −0.5 ≤ a * ≤ 0 and much more evidently from black hole with a * = 0.9375.
The azimuthal profiles derived from the reconstructed images shown in the middle panels of Fig. 5 clearly highlight that the azimuthal variability is observationally detectable.
The comparison between the top and middle panels reveals consistent time development of ϕ obs, peak (t) in both the GRRT and reconstructed images.The time-aggregated distributions of ϕ obs peak (t) are reported in the bottom panels using GRRT (blue histogram) and reconstructed images (red histogram).Note how the negative spin cases exhibit a wide range of peak azimuthal angles, whereas the positive spin cases display sharply peaked histograms.
Figure 6 is similar to Fig. 5, but in it we show the azimuthal variability in the disk region only, r obs < 100 µas (corresponding to r obs /r g < 26).The temporal evolutions of ϕ obs, peak (t) (white points) show more substantial variability than those in the jet region (Fig. 5) and a similar trend is evident in the jet region, where the variations of the white points are larger for the negative spin cases compared to the positive spin cases.
In Fig. 7, we summarize the standard deviation of the azimuthal angle at peak brightness ϕ obs, peak (t) in the jet (circle points) and the disk regions (square points), utilizing both GRRT (blue) and reconstructed images (red).The peak azimuthal angle of the disk region is determined by the gravitational lensing and beaming effects, as well as the brightness distribution of the accretion flow.The transition value be- tween a * ≤ 0 and a * > 0 is approximately σ(ϕ obs, peak ) ∼ 30 • in both regions, which allows us to constrain the direction of the black-hole rotation.The standard deviation of the peak azimuthal angles exhibits similar values across nearly all spin cases in both the jet and the disk images.Moreover, we can see that the disk component σ(ϕ obs,peak ) ∼ 90 • has an azimuthal variation that is larger than that of the jet component (∼ 40 • ) in the case of a * = −0.5, which can be explained with the different rotational direction between the black hole and accretion flow.For a * ≥ 0, the gravitational lensing and beaming effects cause the mean peak position angle to be distributed around ∼ −60 • (see Fig. 6).On the other hand, at a * = −0.9375, the mean peak azimuthal angle is distributed around ∼ 63 • due to the strong frame dragging effect of the counter-clockwise rotation of the black hole, which disrupts the rotation direction of the accretion flow in the opposite direction.In the case of a * = −0.5, the effect of the black hole's frame dragging is not as dominant over the opposing flow of the accretion disk as it is when a * = −0.9375,and thus, large fluctuations occur due to the competition between these two effects.This unique property of azimuthal variation seen in the case of a * = −0.5 has potential utility in distinguishing between observational results for a * = −0.9375,−0.5, and 0.0.

Radial wave propagation of jet inhomogeneities
In this section, we shift our focus to the radial variability in the jet region and propose observational quantities to estimate the magnitude of the black-hole spin.To this scope, we investigate the time development of radial jet propagation.We introduced a distribution of the relative brightness temperature where T B (t, r obs , ϕ obs ) represents the brightness temperature in polar coordinates at the observational time t, and ∆t = 10t g denotes the time interval of the GRRT simulations.We calculated the relative brightness in regions where T B (t, r obs , ϕ obs ) and T B (t + ∆t, r obs , ϕ obs ) exceed the sensitivity threshold of the ngEHT, T B,thres = 2 × 10 7 K.
The GRRT images clearly show the propagation of inhomogeneous wave-like features in the radial direction.In Fig. 8, we report the inhomogeneous radial wave observed in the case of a * = 0.9375.The relative distributions between the neighboring time stamps (top and left three panels) demonstrate the propagation of inhomogeneity in the radial direction.Comparing the relative brightness distributions between the GRRT (top left panels) and the reconstructed images (bottom left panels) provides evidence that similar wave-like features propagating outwards in the radial direction can be detected through expected ngEHT observations.
To better highlight the temporal evolution of the radial wave and measure its speed, we show in the right panels of Fig. 8 the time development of the radial profile at the azimuthal angle of ϕ obs = −20 • .Intermittent wave propagation is clearly observed, as evidenced by the presence of five red curves in the spacetime diagram within the time range of 1200 ≤ t/t g < 1400 based on both the GRRT (top panel) and reconstructed images (bottom panel).Each red/blue stripe exhibits a linear behavior, indicating the intermittent radial propagation of wave inhomogeneity with an almost continuous velocity.This feature sets the stage for our forthcoming examination of the spin dependence on the radial wave velocity in the subsequent section.

Radial velocity of jet inhomogeneity
By utilizing the spacetime diagram in the right panels of Fig. 8, we conducted an investigation into the temporal and radial development of the inhomogeneous waves in the jet region.We estimated the time and azimuthal distribution of the number of waves N , relative brightness f , and radial velocity v r of each wave through the following steps: 1. Given an azimuthal angle ϕ, we collected radial waves with f (t, r obs ) > 0.1, that is, with an excess of 10% in the relative temperature brightness.Continuity of In this table, σjet(ϕ obs ) and σ disk (ϕ obs ) denote the standard deviation of the peak azimuthal angles in the jet and the disk regions, respectively (refer to Fig. 7).Furthermore, v r true represents the maximum radial wave velocity along the azimuth (refer to Fig. 11).The values enclosed in parentheses are the estimates derived from the reconstructed images.wave propagation is defined as a set of temporally and radially continuous points, spanning a time duration of more than 30 t g .Based on this process, we obtained a set of events (t, r obs ) and the corresponding relative brightness f of each wave train.
2. We repeated this procedure at each azimuthal angle, resulting in a set of waves (t, r obs , ϕ obs , f ) j , where i denotes the wave index 1 ≤ j ≤ N tot and N tot is the total number of detected waves.
3. We conducted a linear fitting to each wave j using the least square method, such that the slope of the linear function r obs, fit (t j ) provides the radial velocity v r .Note that because the estimated velocity is projected to the image plane, the deprojected relativistic velocity v r true is calculated with the observed apparent one v r (e.g., Falla & Floyd 2002): 4. Finally, the velocity was averaged both temporally and azimuthally across sections with widths ∆t and ∆ϕ.
For the formers, we selected ∆t = 200t g , which corresponds to the propagation time in the radial jet region (100 µas ≤ r obs < 300 µas) when v r ≥ 0.1.
For the latter, we have chosen an azimuthal section, ∆ϕ = 30 • , which encompasses the range spanned by azimuthal jet brightness for a * ≥ 0.5 (see Fig. 7), and is consistent with the width of the bottom histograms in Fig. 5.The validity of the radial velocity estimate was inspected using geometric jet models shown in Appendix A.
Following steps (1)-(3), we investigated the radial wave properties in the jet region, and Fig. 9 shows the number of detected waves N with respect to the relative brightness f .The top and middle panels show the distributions of N/max(N ) at the residual radius r obs − r obs, fit and the relative brightness of each wave, utilizing the GRRT (top panels) and reconstructed images (middle panels), respectively.Clearly the maximum value of the relative brightness f depends on the absolute value of the spin, and for N/ max(N ) = 0.1 (white curves), f extends to approximately 1.5, 1.4, 1.4, 1.3, and 1.6 in the cases of a * = −0.9375,−0.5, 0.0, 0.5, and 0.9375, respectively.This behavior suggests the inhomogeneity of the radial waves increases with the magnitude of the black-hole spin.The number of radial waves, N tot (left value), monotonically increases with the absolute value of the spin in the case of the GRRT images (top panels), and a similar trend is seen also for the reconstructed images (middle panels).However, they do not quite accurately capture the spin dependence observed in the GRRT images due to the uncertainties in the synthetic observations.The bottom panels show the wave width of the fit lines.Each component is radially symmetric and has a FWHM of ∼ 20 µas, which is the restoring beam size of the image (see Section 2.2).The symmetric profiles, characterized by the FWHM, indicate the accurate detection of the center of jet waves using linear fitting.This also provides insight into the quality of the fitting process, with the validation detailed in Appendix A.
Figure 10 presents the estimated radial velocities v r true for each spin.The top and middle panels display the distributions of v r true at each azimuthal angle and time using the GRRT and reconstructed images, respectively.The estimated velocities exhibit significant variations in both time and azimuthal directions for a * < 0.5, while demonstrating more robust statistics for a * ≥ 0.5.This behavior is consistent with the spin dependence of the azimuthal variability of the jet discussed in Section 4.1.The bottom panels show the median velocity with respect to the temporal direction for each azimuthal angle, where the error bars correspond to the median absolute deviations.Note that the velocity peaks around the rotation axis of the black hole (ϕ = 18 • ) for all spins.Finally, Fig. 11 summarizes many of the results presented so far and the relationship between spin and the maximum radial velocity in the azimuthal direction.The estimated velocity, as indicated by the GRRT (blue) and reconstructed (red) images, demonstrates a monotonic dependency on the absolute spin value.These estimates allow us to distinguish between the cases of a * = 0.5 and a * ≥ 0.9375.As discussed in Section 4.1, the azimuthal variability of the jet and the disk enables differentiation in the case of a * ≤ 0 (the orange region in Fig. 11).As a result, Fig. 11 provides evidence that the azimuthal and radial variations of the jet and the disk depend on the black-hole spin, thus enabling us to extract spacetime information (Table 4).Similar spin dependences are also observable at different inclination angles i, with the measurement of radial velocity and azimuthal variation providing constraints on a * and i (Appendix B). .The number of detected waves with respect to the relative brightness f across the entire time and azimuthal domains (from left to right: a * = −0.9375,−0.5, 0.0, 0.5, and 0.9375).The top and middle panels show the distributions of N/max(N ) at the residual radius r obs − r obs, fit and the relative brightness of each wave, utilizing the GRRT and reconstructed images, respectively.Here, r obs, fit represents the fitting radius of each wave in the spacetime diagram, the white curve represents the boundary where N/max(N ) = 0.1, and the leftmost number Ntot denotes the total number of detected waves.The bottom panels depict the vertical integration of N from the upper panels, with normalization to its maximum value.

PHYSICAL INTERPRETATION OF JET RADIAL WAVES
The origin of radial jet waves can primarily be attributed to the motion of relativistic plasma along the jet, a phenomenon both confined and accelerated by the magnetic field.By examining the fluid velocities (|v| = √ v i v i , where v i and v i are the covariant and contravariant three velocities of the accretion flow) and Alfvén velocities (v Alf = σ B /(σ B + h)), we investigate the origin of the radial jet wave velocity derived in the previous section.Here, σ B = b 2 /ρ is the plasma magnetization, b is the norm of the magnetic field in the fluid frame, and ρ is the rest-mass density.Figure 12 reports radial distribution of the velocities averaged over time and around the jet; the red and blue curves denote the fluid and Alfvén velocities, respectively.The dashed and solid lines indicate the average values around the jet axis (σ B > 3) and in the jet region (σ B < 3 and Be > 1.02), respectively.In each region and each spin case, the fluid velocity (blue curves) exhibits a monotonic increase with the radius due to the diminishing gravitational attraction.Conversely, the Alfvén velocity (red curves) decreases monotonically with the radius as it moves farther from the center where the magnetic field is most concentrated.
Both velocities have a monotonic dependency on the spin around the jet axis (see dashed curves).The fluid and Alfvén velocity in the region averaged in the radial direction, are given by |v| = 0.81, 0.74, 0.70, 0.82, 0.87 and v Alf = 0.85, 0.93, 0.94, 0.93, 0.87 , corresponding to the blackhole spins of a * = −0.9375,− 0.5, 0.0, 0.5, and 0.9375, respectively.These results indicate that the accretion flows around the jet region accelerated with the strong magnetic fields represented by the Alfvén velocity (v Alf ∼ 0.90) and the fluid velocity monotonically increases with the magnitude of the spin.
We further investigate the spin dependency of the fluid and Alfvén velocities within the jet region.In Fig. 13, we present a comparison between the spatially averaged fluid and Alfvén velocities in the jet region (red and blue points) and the estimated velocity derived from the jet radial waves (Section 4.3, depicted in black).The accompanying error bars for the fluid and Alfvén velocities represent the standard deviation.The fluid velocities align closely with the estimates from the jet radial waves.In contrast, the Alfvén velocities do not exhibit a clear spin dependency.It is further noted that the fluid velocity in the jet region is predominantly radial, with ratios of the radial three-velocity √ v r v r and |v| averaged in the jet region are 85 %, 90 %, 93 %, 92 %, 95 % for a * = −0.9375,−0.5, 0.0, 0.5, and 0.9375, respectively.The results indicate that inhomogeneity within the jet serves as an effective marker, tracing the contribution of spin to the relativistic fluid motion.

SUMMARY AND FUTURE PROSPECTS
We investigated the horizon-scale dynamics of the jet and the accretion disk and demonstrated the potential of blackhole spin measurement.We focused on the best-bet GRMHD models whose continuum spectra are consistent with that of the M87 * observations within a small set of model parameters (Fromm et al. 2022;Cruz-Osorio et al. 2022).The azimuthal variation of the jet and disk brightness reveals information about the spin direction, with the azimuthal angle of peak brightness showing temporal stability in cases of a * > 0 and instability when a * ≤ 0. Furthermore, this variation has the potential to distinguish negative spin cases.In the jet region (100 µas ≲ r obs ≲ 300 µas), inhomogeneous waves propagate radially with relativistic velocities (0.26 ≲ v r true ≲ 0.50).The radial wave velocity, as estimated from the GRRT images, shows a monotonic dependence on the absolute value of the black-hole spin.The reconstructed images from the ngEHT have demonstrated the detectability of these dynamics and the accuracy of the estimated black-hole spin.This paper provides the essence of the spacetime measurement with the horizon-scale dynamics and future theoretical and observational projects.
A survey with a broader library of GRMHD simulations will provide a more comprehensive view of spin measurement under various possibilities of spacetime, magnetic field, .Spin dependency of the azimuthally maximum value of radial wave velocity v r true (refer to the bottom panels in Fig. 10), where the blue and red points represent the estimated values obtained from the GRRT and reconstructed images, respectively.The top labels in the orange and cyan-colored regions indicate potential characteristics that can be used to constrain the spin within each spin range.For clarity, the red points are slightly shifted to the right (0.05) in the horizontal direction.accretion flow properties, and inclination angles.For a survey with an expanded library of GRMHD simulations, it is essential to investigate each regime of parameter constraints with individual properties: continuum spectrum, visibility fitting, polarization, and jet-disk dynamics (Fromm et al. 2022;Event Horizon Telescope Collaboration et al. 2019f, 2021a;Davelaar et al. 2023).In addition to the expansion of model parameters, the dynamics survey has the potential to provide additional information on tilted disk parameters depending on the Lense-Thirring precession of the accretion disk and jet (e.g., Liska et al. 2018;Chatterjee et al. 2020).We also focused on the simplest case without gain uncertainty and performed image reconstruction.To investigate the fiducial morphology under more realistic conditions, we need to conduct imaging surveys with numerous imaging parameters (Event Horizon Telescope Collaboration et al. 2019d, 2022c) with the more realistic data including gain uncertainties, polarization leakages, and instrumental effects.
Alongside the prospective utility of ngEHT observations, it is worth noting that the annual azimuthal variation in the photon ring, as observed through current and impending EHT studies, is also synergetic with our proposed method.The annual EHT observations, which offer less dense (u, v) cov-erages compared to ngEHT, are either under analysis or in the planning stages to investigate variations of the photon ring.Particularly, RML imaging, as demonstrated with the 2017 EHT data, has furnished an azimuthal angle assessment of the photon ring, achieving an accuracy of approximately ∼ 20 • (Event Horizon Telescope Collaboration et al. 2019d).By continuous observations of the azimuthal angle variabilities and capturing the temporal evolution of images within the disk region, discerning differences in the azimuthal angle variation correlated to the black hole's spin direction (roughly 60 • ) will become attainable.Conversely, the task of differentiating small differences in azimuthal angle variabilities between (a * = 0.5) and (a * = 0.9375) (around 2 − 3 • , as outlined in Table 4) is challenging.The observation of the jet's temporal progression, augmented by the denser (u, v) coverage facilitated by ngEHT observations, will be an essential requirement for constraining the magnitudes of the spin.
The quality of horizon-scale images has improved significantly due to recent advancements in imaging techniques (e.g., Chael et al. 2022;Mueller & Lobanov 2022;Arras et al. 2022).These developments potentially pave the way for the application of our spin measurement methodology to EHT observations from 2021 onwards.These results, in conjunction with EHT observations, are expected to synergize with current and future multi-wavelength observations (Lu et al. 2023).By utilizing our methodology in future EHT observations, theoretical and observational studies will be able to achieve a much more comprehensive measurement of the black-hole spacetime of M87 * .(https://achael.github.io/eht-imaging/)(Chael et al. 2018b), and SMILI (https://github.com/astrosmili/smili)(Akiyama et al. 2017a,b).

APPENDIX A. INSPECTION OF RADIAL VELOCITY ESTIMATES WITH GEOMETRICAL JET MODELS
In Section 4.3, we introduced a method for estimating the radial wave velocity of jet inhomogeneity.To validate these velocity estimates, we applied the method to simple geometric jet models and evaluated the accuracy of the resulting velocity estimates.The intensity distribution of the geometric model is given by I(t, r, ϕ) = I 0 (t)G(t, r, ϕ) exp − r − v r t ∆r deviation of one.The intensity I 0 (t) is defined such that the total flux equals 1.0 Jy at each time.
Figure 14 compares the estimated velocities to the input (true) velocities, showcasing typical examples with ∆r = 20 µas and ∆ϕ = 40 • .The velocity estimates are derived from the input images following the same strategy described in Section 4.3.Across the input velocity space in this paper (0.2 ≤ v r ≤ 1), the estimated velocities accurately reproduce the ground-truth values with an error of about ∆v r ∼ 0.01.

B. INCLINATION ANGLE DEPENDENCY
The main text focuses on the best-bet model parameters, which reproduce the continuum spectrum of M87 * at a fixed viewing angle (i = 160 • ).In this appendix, we demonstrate that the broad spin dependency of both the azimuthal variability of the jet (Section 4.1) σ(ϕ obs, peak ) and radial wave velocity v r true (Section 4.3) is observable at different inclination angles.
In Fig. 15, we show the inclination angle dependency of the azimuthal variation σ(ϕ obs, peak ) and the radial wave velocity v r true for each spin (the red, blue, and green curves represent the cases of |a * | = 0.9375, 0.5, and 0, respectively) with the GRRT images.The transition value of σ(ϕ obs, peak ) between a * ≤ 0 and a * > 0 is σ(ϕ obs, peak ) ∼ 30 nation angles.Thus, the measurement of the radial velocity and azimuthal variability could potentially constrain both the spin value and the viewing angle.

Table 1 .Figure 1 .
Figure1.Aggregate baseline coverage for the ngEHT observation with the expected telescopes(Roelofs et al. 2023).The blue and red data points represent the (u, v) coverage obtained from the telescopes integrated in 2017 and the ngEHT, respectively, with corresponding observational time durations of 7 hours and 24 hours.The axes of (u, v) plane are in units of giga-wavelengths at the spatial frequency of 230 GHz.

Figure 2 .
Figure 2. Comparison between the top panels displaying snapshot images obtained from the GRRT simulation and the middle panels demonstrating reconstructed images using synthetic data based on the expected ngEHT observations at a frequency of 230 GHz (Section 4), specifically focusing on the case characterized by a * = 0.9375.The images in the top right and middle right panels are restored using a circular Gaussian beam with a full width at half maximum (FWHM) of 20 µas, as indicated in the lower right corners.The unit represented in each image is the logarithm of the brightness temperature, denoted as TB.The bottom panels depict the ratio between the GRRT and reconstructed images, with the bottom left panel displaying the ratio without restoration, and the bottom right panel representing the ratio after restoration with the Gaussian beam.
Figure3.The time-averaged properties of the accretion disk and the jet at 230 GHz.The left panel displays the time-averaged GRRT image in spherical polar coordinates for the case of a spin and inclination angle of (a * , i) = (0.9375, 160 • ).The green and cyan curves represent the disk and jet+disk regions, respectively, which can be detected with the typical ngEHT sensitivity (= 2 × 10 7 K).In the middle panel, the radial profile of the azimuthal maximum brightness temperature of the jet+disk is shown.The expected sensitivity of the 2017 EHT and ngEHT observations is presented as horizontal gray lines.In the right panel, the same format is followed as in the middle panel, but for the disk emission defined by the Bernoulli parameter (Be < 1.02).The vertical gray lines in the middle and right panels indicate a radius of 100 µas, which is the typical radius of the disk region for each spin case.

Figure 4 .
Figure4.The spin dependency of the 230 GHz light-curves and power spectral densities (PSDs).From left to right, the different panels refer to black holes with a * = −0.9375 to 0.9375.Top panels: the light curves emitted from both jet and disk (black), jet-only (blue), and disk-only (red) regions.Bottom panels: the power spectral densities of the light curves from the upper panels, emitted from each region, where τ = 2000tg is the entire simulation time window.in line with the analysis in Appendix B of Event Horizon TelescopeCollaboration et al. (2019e).The ratios of the total flux density are 66 %, 85 %, 79 %, 87 %, and 81 % for a * = −0.9375,−0.5, 0.0, 0.5, and 0.9375, respectively, whereas the MAD thermal models yield a value of about 90 % (see Appendix B in Event Horizon TelescopeCollaboration et al. 2019e).This discrepancy can be attributed to the impact of non-thermal particles, which result in a comparable continuum spectrum spanning from the radio to the near-infrared range and influence the jet morphology at a frequency of 86 GHz in the case of M87 * (seeFromm et al. 2022;Cruz-Osorio et al. 2022).The variability characteristics of each light curve are also apparent in the power spectral densities (PSDs) P ω , where ω denotes the angular frequency (bottom panels of Fig.4).The frequency profiles of the PSDs, that are characterized by P ω ∝ ω −2.4±0.5 , are broadly consistent with those produced by red noise (P ω ∝ ω −2 ).This property was already reported in previous research that utilized a large library of GRMHD simulation models without non-thermal particles(Georgiev et al. 2022).Interestingly, some PSDs display a feature that could be associated with a quasi-periodic oscillation (QPO).This trend is especially pronounced for a * = 0.9375, where the PSD exhibits peak components at time scales of 1/ω ∼ 320t g and

Figure 6 .Figure 7 .
Figure6.The same as in Fig.5, but each profile was radially averaged only in the disk region r obs < 100 µas.

Figure 8 .
Figure 8.The radial propagation of wave inhomogeneity within the jet region in the case of a * = 0.9375.The top and bottom panels represent the distributions using the GRRT and reconstructed images, respectively.The three panels on the left-hand side show the polar distribution of temperature inhomogeneity at corresponding timestamps (t/tg = 1280, 1290, and 1300).The panels on the right-hand side display spacetime diagrams illustrating the inhomogeneity distributions at an azimuthal angle of ϕ obs = −20 • .
Figure9.The number of detected waves with respect to the relative brightness f across the entire time and azimuthal domains (from left to right: a * = −0.9375,−0.5, 0.0, 0.5, and 0.9375).The top and middle panels show the distributions of N/max(N ) at the residual radius r obs − r obs, fit and the relative brightness of each wave, utilizing the GRRT and reconstructed images, respectively.Here, r obs, fit represents the fitting radius of each wave in the spacetime diagram, the white curve represents the boundary where N/max(N ) = 0.1, and the leftmost number Ntot denotes the total number of detected waves.The bottom panels depict the vertical integration of N from the upper panels, with normalization to its maximum value.

Figure 10 .
Figure10.Azimuthal and temporal distributions of the radial wave velocities v r true for each spin.The top and middle panels show the distribution of radial velocity using the GRRT and reconstructed images, respectively.The lower panels present the median values of the radial velocity in the temporal direction, derived from both the GRRT (blue) and reconstructed (red) images, with the error bars indicating the corresponding median absolute deviations.
Figure11.Spin dependency of the azimuthally maximum value of radial wave velocity v r true (refer to the bottom panels in Fig.10), where the blue and red points represent the estimated values obtained from the GRRT and reconstructed images, respectively.The top labels in the orange and cyan-colored regions indicate potential characteristics that can be used to constrain the spin within each spin range.For clarity, the red points are slightly shifted to the right (0.05) in the horizontal direction.

Figure 12 .Figure 13 .
Figure 12.Radial distribution of the fluid velocity (|v|, red curve) and Alfvén velocity (v Alf , blue curve) averaged over time and around the jet for each spin value.The dashed and solid curves correspond to the average around the vertical jet axis (σB > 3) and in the jet region (Be > 1.02 and σB < 3).

2(Figure 14 .
Figure14.The radial velocity of the geometric jets as defined in Section A. The gray curve represents the input (or 'true') radial velocity, while the red curve represents the estimated radial velocity.

Figure 15 .
Figure 15.The inclination angle dependence of the estimated azimuthal variation (top panel) and radial velocity (bottom panel).For clarity, the curves in the case of a * ≥ 0 are slightly (1 • ) shifted in the horizontal direction.

Table 4 .
Summary of the spin dependency of the disk-jet dynamic.