M 87: a cosmic laboratory for deciphering black hole accretion and jet formation

A fraction of galaxies emit an immense amount of energy from the compact central regions, often surpassing the combined brightness of their entire host galaxies. These are known as active galactic nuclei (AGNs), and it is widely accepted that they are powered by the accretion of material onto supermassive black holes (SMBHs) at their cores (Lynden-Bell 1969; Rees 1984). The quest to understand the nature of AGNs and SMBHs has persisted since their initial discovery and remains a paramount focus of the current astrophysics. Furthermore, there is growing evidence that the energetic activities of AGN profoundly influence the formation and cosmological evolution of galaxies and galaxy clusters, as represented by the empirical correlation between the stellar velocity dispersion of galaxy bulges and the mass of central SMBHs (Magorrian et al 1998; Kormendy and Ho 2013). This feedback mechanism may manifest through the intense radiation emitted by luminous accretion disks or the mechanical energy of powerful jets/outflows emanating from the nuclei.

AGNs display diverse observational properties, with various types and subclasses distinguished based on different features. While some of this diversity can be attributed to an orientation effect, where the central engine remains essentially identical, intrinsic differences in the central object also contribute to the range of AGN activity. Most notably, AGNs are divided into two distinct classes based on the ratio of radio to optical luminosity (i.e., radio loudness): radio-quiet (RQ) AGNs and radio-loud (RL) AGNs (Antonucci 1993; Urry and Padovani 1995). The former constitutes the predominant population of AGNs, with their spectral energy distributions (SED) primarily dominated by thermal emission from accretion disks observed across optical, ultraviolet, and X-ray wavelengths. In contrast, the latter class represents only 10–15% of AGNs (Kellermann et al 1989), but is characterized by powerful jets consisting of relativistically beamed outflows composed of highly magnetized plasma. The relativistic jets make the appearance of AGNs more spectacular. When viewed from large orientation angles, radio jets extend from kiloparsec (kpc) to Megaparsec (Mpc) scales, displaying a two-sided morphology centered on the nucleus. Such misaligned RLAGNs are referred to as radio galaxies. Conversely, when the jets are viewed from small angles, the apparent structure of the source becomes compact but shows intense nonthermal radiation and variability from radio to $\gamma$-rays due to a strong Doppler-boosting effect. This type of jetted AGNs is known as blazars. Both radio galaxies and blazars are further classified into subclasses based on the observed jet powers, radio morphology and the characteristics of the observed SED and emission lines (Urry and Padovani 1995).

Despite the wealth of observational features seen across the entire electromagnetic spectrum, our understanding on the physical processes governing the accretion and ejection associated with SMBHs is still incomplete: How exactly is matter transferred from galactic scales to event horizon scales? What control the power of relativistic jets? What causes the intrinsic diversity of AGNs? One obvious challenge is that AGN and its environment are highly complex and a variety of complicated physical processes are at work. Besides, the relevant physical scales associated with the BH accretion/ejection are tiny, typically within parsec (pc) scales or much less, making them difficult to spatially resolve with most current astronomical facilities. Hence, the key to advancing our understanding lies in identifying a nearby target that enable detailed spatial resolution and investigation of the SMBH vicinity, coupled with the availability of high-resolution observational facilities.

In this review we specifically focus on M 87 (NGC 4486), a nearby giant elliptical galaxy in the center of the Virgo cluster located at a distance of $D=16.8$ Mpc (Blakeslee et al 2009; Bird et al 2010; Cantiello et al 2018). This object was originally discovered and cataloged by Charles Messier in 1781. In 1918, Heber Curtis at Lick Observatory detected ‘a curious straight ray …apparently connected with the nucleus by a thin line of matter’ in their optical image of this galaxy (Curtis 1918). Although the exact physical origin of this feature was not known at that time, this was later recognized as the first discovery of an astrophysical jet. About 40 years later, Walter Baade discovered that the light from the ‘jet’ was strongly polarized, first indicative of nonthermal synchrotron origin of the jet emission (Baade 1956). M 87 has also been known as one of the brightest radio sources on the sky (as Virgo A or 3C 274) since the early days of radio astronomy (e.g., Bolton et al 1949). The large-scale radio structure associated with the galaxy extends as large as 80 kpc (Mills 1952; Baade and Minkowski 1954; Owen et al 2000), and Fanaroff and Riley classified this source as ‘Class I’ based on its radio morphology and relatively low radio luminosity (Fanaroff and Riley 1974)¹¹1Based on more than 50 samples of radio galaxies, Fanaroff and Riley (1974) classified radio galaxies into Class I and Class II, which are currently known as FR-I and FR-II, respectively. FR-I radio galaxies often have symmetric radio jets whose intensity falls away from the nucleus, while FR-II radio galaxies typically exhibit more highly collimated jets leading to bright hot spots far from the nucleus. Fanaroff and Riley (1974) also found that FR-I/FR-II sources have lower/higher luminosities, with a division around $\sim 2\times 10^{25}\,{\rm WHz^{-1}str^{-1}}$ at 178 MHz. These facts have led to the suggestion that FR-I sources host relatively weakly accreting AGNs, while FR-II sources contain more powerful AGNs.. After 1980s the advent of the Very Large Array (VLA) revealed the kpc-to-subkpc-scale radio structures in great detail (along with a well-collimated one-sided jet) thanks to its high-dynamic-range imaging capability at subarcsecond-scale angular resolution (e.g., Biretta et al 1983; Owen et al 1989). In 1990s, high-resolution optical spectroscopic observations with the Hubble Space Telescope (HST) revealed a rotating ionized gas disk in the central pc-scale nuclear regions of the galaxy (Harms et al 1994; Ford et al 1994; Macchetto et al 1997), which provided the first strong hint for the presence of a central black hole with exceeding a billion solar masses in the center of this galaxy (see Sect. 3 for details).

Since its discovery, M 87 has been a privileged target for studying the physics of accretion and jet formation associated with an SMBH. Compared to distant quasars, M 87 is much closer to us. Combined with a large mass of the central SMBH, this makes the apparent diameter of the SMBH second largest on the sky after the Galactic Center Sgr A*, which allows us to resolve the source structure into smaller gravitational scales at a given angular resolution. On the other hand, unlike Sgr A*, M 87 exhibits a powerful relativistic jet, which is the most outstanding feature characterizing RLAGNs. Moreover, M 87 is sufficiently bright across the entire electromagnetic spectrum not only in radio/optical but also in X-rays and $\gamma$-rays, making this source an excellent target to study broadband multi-wavelength (MWL) properties associated with AGN activities (Fig. 1).

In particular, rapid and remarkable advances have recently been made by high-resolution very-long-baseline-interferometry (VLBI) observations at radio wavelengths. The most groundbreaking example is the detection of the black hole shadow with the Event Horizon Telescope (EHT; Event Horizon Telescope Collaboration et al 2019a, b, c, d, e, f), a global VLBI network operated at 1.3 millimeters (mm). This provided the first visual evidence for the existence of SMBH as the central engine of AGN (see Sect. 3). Global VLBI observations at mm wavelengths (including 3 mm) are also capable of probing the inner part of accretion flows for M 87 (see Sect. 4). At even longer wavelengths, VLBI observations are generally more sensitive to the larger-scale emission, enabling us to image the collimation and acceleration regions of relativistic jets extending from near-horizon scales to galactic/kpc scales (see Sect. 5).

In summary, M 87 provides a rare opportunity to address various key questions relevant to (RL)AGNs as well as the physical connection among the SMBH, accretion flows and relativistic jets, especially under the low accretion/Eddington regime (see Sect. 4), across all electromagnetic wavelengths and spatial scales. In fact, despite being a single object, M 87 has been studied for many decades in the astrophysical community, with international workshops focused on this object (e.g., Blandford 1999). Attention to this object has been growing even more rapidly in the last few years. Fig. 2 displays annual trends in article publication related to M 87. While the number of annual publications has steadily increased from 1970s to 2010s, one can see a dramatic rise in publications since 2019/2020, very likely motivated by the first EHT imaging results. Therefore, now it would be a good time to review the current observational understanding of this active galaxy.

In this review, we provide a comprehensive overview of the current observational understanding of the active galaxy M 87, covering a broad range of spatial scales from event horizon scales to kpc scales. The review is structured as follows. In Sect. 2, we first summarize the key parameters characterizing AGN/SMBH intrinsic activity. In Sect. 3, 4 and 5, we highlight various observational progress on M 87, dividing the topics into three key components: the central SMBH, the accretion flows onto SMBH, and the relativistic jet, respectively. In Sect. 6, we summarize our current understanding of this source and highlight some future prospects. Finally, we conclude the review in Sect. 7. Readers may also refer to Boccardi et al (2017) for a focused review on millimeter-VLBI observations of AGNs, and to Blandford et al (2019) for an extensive review covering both observational and theoretical efforts on studies of AGN jets in general.

SMBHs are found in most, if not all, galactic nuclei with a mass range of $M_{\rm BH}$, from $\sim$$10^{6}\,M_{\odot}$ to $\sim$$10^{10}\,M_{\odot}$ (Kormendy and Ho 2013). RLAGNs exist in a variety of host galaxies, from disk-dominated spirals to giant ellipticals, with a similar $M_{\rm BH}$ range as described above (Ho 2002). It seems that there is no clear difference in the BH mass range ($10^{6}-10^{10}\,M_{\odot}$) between RLAGNs and RQAGNs (Woo and Urry 2002), implying that $M_{\rm BH}$ may not be a dominant factor in characterizing AGN activity at radio bands.

Blazars, the other population of RLAGNs, fall into two sub-classes i) flat-spectrum radio quasars (FSRQs; objects with emission-line dominated spectra) and ii) BL Lac objects (BLOs; nearly line-less objects), with generally higher SED peaks (for both synchrotron and inverse Compton radiation) in FSRQs than in BLOs. The $M_{\rm BH}$ of blazars spans from $\sim$$10^{7.5}\,M_{\odot}$ to $\sim$$10^{9.5}\,M_{\odot}$ (Wu et al 2002; Shaw et al 2012; Xiao et al 2022) with a possible localization for BLOs at $10^{8}-10^{9}\,M_{\odot}$ and for FSRQs at $10^{8}-10^{9.5}\,M_{\odot}$ (Ghisellini et al 2010; Chen et al 2015).

The majority of RLAGNs (‘jetted’) have been believed to be hosted by elliptical galaxies with high BH masses $M_{\rm BH}\gtrsim 10^{8}\,M_{\odot}$ (Fabian and Canizares 1988), while late-type galaxies, such as spirals with lower BH masses $M_{\rm BH}\lesssim 10^{8}\,M_{\odot}$, were traditionally considered as RQAGNs (‘non-jetted’) (Laor 2000). This conventional hypothesis, however, has been updated during the last two decades by discovering radio-loud narrow-line Seyfert 1 galaxies (RLNLSy1s) with $M_{\rm BH}\sim 10^{6}-10^{8}\,M_{\odot}$ (e.g., Komossa et al 2006; Yuan et al 2008; Mathur et al 2012). Furthermore, the detection of $\gamma$-ray emitting RLNLSy1 objects (Abdo et al 2009) has opened a new era for understanding the jetted RLAGNs hosted in diverse types of galaxies, including spirals, with a wide range of $M_{\rm BH}\sim 10^{6}-10^{10}\,M_{\odot}$.

There are various independent methods proposed to calculate $M_{\rm BH}$: the traditional virial BH mass, the stellar velocity dispersion or the bulge luminosity, and the variation time scale. The $M_{\rm BH}$ of M 87 has been estimated with several methods as introduced in Sect. 3.

BHs can rotate like other stars. Their rotation speeds could be close to the speed of light, causing a dramatic change in the metric and in the way matter moves (the mass accretion and the energy extraction). The properties of a BH are uniquely characterized by its mass, angular momentum, and electrical charge, based on the no-hair theorem of GR. In reality, a charged BH may be hard to exist due to its neutralization by discharge to its surroundings (difficult to maintain particle creation and/or accretion of matter onto the hole). Therefore, astrophysical BHs can be determined solely by the mass and angular momentum (henceforth referred to as spin).

An uncharged, rotating (spinning) BH is described by the Kerr metric (Kerr 1963) with the BH mass $M_{\rm BH}$ and angular momentum $J$. A dimensionless angular momentum (spin) parameter ($-1\leq a\leq 1$) that measures a fraction of the angular momentum $J$ relative to its maximum value is defined as follows:

where $G$ and $c$ are the gravitational constant and the speed of light, respectively. Here, we introduce the ‘gravitational radius’ $r_{\rm g}\equiv GM_{\rm BH}/c^{2}$ as the standard definition. The BH can be rotating in a prograde ($a>0$) or retrograde ($a<0$) manner in that system.

We also do want to locate the event horizon ($r_{\rm H}$), the boundary of the BH. In the Schwarzschild metric (a static metric for the spacetime around a non-rotating BH with $a=0$, where the spatial coordinates are now spherical-polar, Schwarzschild 1916), we have $r_{\rm H}=2GM_{\rm BH}/c^{2}$ ($r_{\rm S}$: Schwarzschild radius). But for the Kerr metric, the radius of the BH event horizon is described as $r_{\rm H}\equiv r_{\rm g}(1+\sqrt{1-a^{2}})$, where $r_{\rm H}\rightarrow r_{\rm g}$ with maximally spinning cases in both prograde and retrograde fashions. Note that the horizon is still spherical even for a spinning BH. Another important radius appears in the Kerr metric, $r_{\rm E}\equiv r_{\rm g}(1+\sqrt{1-a^{2}\cos^{2}{\theta}})$, which is larger than $r_{\rm H}$ for all $\theta$ and $\phi$ (azimuthal and polar angles, respectively), except at the poles where the two are equal. Note that at the equator, always $r_{\rm E}=2r_{\rm g}$, for any BH spin.

The region $r_{\rm H}<r<r_{\rm E}$, called the ‘ergosphere’, is peculiar: everything must rotate in the same direction as the BH, just as everything must fall toward the BH within the horizon. The BH rotation causes the inertial frames to be dragged at an angular frequency $\Omega\approx\Omega_{\rm H}/(r/r_{\rm H})^{-3}$, where

is the angular frequency of the event horizon Misner et al (1973). This implies that the rotation rate $\Omega$ does not vanish beyond the ergosphere, but rapidly drops off outward. Thus, the frame dragging is almost concentrated near or within this region.

In the vicinity of a spinning BH, magnetic field lines advected with the accreting matter, can also rotate. When the field lines thread the event horizon, significant amounts of energy and angular momentum are extracted from the rapidly spinning BH along all field lines (see Fig. 3). Each field line is wound in such a way that the toroidal (azimuthal) field component is generated. An angular frequency of the field line $\Omega_{\rm F}$ will take an intermediate value between $\Omega=\Omega_{\rm H}$ and $\Omega\sim 0$ outside the ergosphere. In general, BHs are not considered to be a perfect conductor so that the field lines threading the event horizon can slip with respect to the horizon: $\Omega_{\rm F}<\Omega_{\rm H}$. Typically, $\Omega_{\rm F}\simeq 0.5\Omega_{\rm H}$ is considered Blandford and Znajek (1977); MacDonald and Thorne (1982), but it depends on the geometry of the BH magnetosphere and/or the spin parameter $a$ Komissarov (2001); Penna et al (2013); Penna (2015).

The mass accretion rate $\dot{M}\equiv dM/dt$, i.e. the rate at which the mass of the accreting compact object increases with time, is a fundamental parameter, and we evaluate it in units of $M_{\odot}\,{\rm yr^{-1}}$ in the astrophysical context. On galactic scales, the hot X-ray emitting intergalactic medium (IGM) through quasi-spherical accretion (Bondi 1952) is widely used to estimate the amount of the accretion power by the SMBH. The Bondi accretion radius $r_{\rm B}$ where the gravitational potential of the SMBH dominates over the thermal energy of the IGM gas, is given by

where $c_{\rm s}$ and $k_{\rm B}$ are the adiabatic sound speed of the gas at the accretion radius (Allen et al 2006) and the Boltzmann constant, respectively. Typical values of the IGM temperature $T$ and $M_{\rm BH}$ are adopted. Correspondingly, the Bondi accretion rate $\dot{M}_{\rm B}$ is given by

for an adiabatic index $\gamma=5/3$ and the mass density $\rho$ with the electron number density $n_{\rm e}$ (Russell et al 2013).

Suppose that the BH accretion starts around the sphere of influence ($\simeq r_{\rm B}$), where the BH’s gravity becomes dominant relative to that of the host galaxy, at an accretion rate similar to the Bondi rate $\simeq\epsilon_{\rm acc}\dot{M}_{\rm B}$ (where $\epsilon_{\rm acc}\leq 1$ is the accretion efficiency). For the Bondi accretion of M 87, readers can refer to Sect. 4.1. In the radiatively inefficient regime, we assume a typical value of the viscosity parameter $\alpha\sim 0.1-0.3$ in the hot, advection-dominated accretion flows (ADAF: $\rho(r)\propto r^{-3/2}$, Narayan and Yi 1995; Mahadevan 1997) and we can refer to the solution of a ‘giant’ ADAF which provides $\epsilon_{\rm acc}\gtrsim 0.5$ (Narayan and Fabian 2011). A flatter radial mass density profile (than the ADAF solution) $\rho(r)\propto r^{-3/2+s}\,(0<s\leq 1)$ suggests a reduction in the mass accretion rate towards the BH: $\dot{M}(r)\propto r^{s}$ (while $\dot{M}(r)=const.$ in ADAF). It would be worth considering that the hot gas is convectively unstable (Abramowicz et al 2002), moderately magnetized, and/or the outflow (wind) co-exists (Blandford and Begelman 1999), as examined by both hydro and magnetohydrodynamic (MHD) simulations (e.g., Yuan and Narayan 2014, references therein).

An extensive numerical survey of slowly rotating magnetized BH accretion through the Bondi radius provides a power-law slope of $s\simeq 0.2-0.5$ Pang et al (2011). Indeed, inside the Bondi radius, a self-consistent profile of the electron number density is obtained in X-ray observations toward M 87 (see Sect. 4). Three-dimensional general relativistic magnetohydrodynamics (GRMHD) simulations provide $\dot{M}\propto r^{0.5}$ at $r\lesssim 100\,r_{\rm S}$ (Sądowski et al 2013; Yuan et al 2015). Considering theoretical/numerical and direct observational results, a substantial reduction of the mass accretion rate compared with the Bondi rate would be expected towards the inner scales.

In the radiatively inefficient accretion flows (RIAFs), the mass accretion rate relative to the Eddingtion rate, $\dot{m}=\dot{M}/\dot{M}_{\rm Edd}$ is well below the unity, where

with the standard radiative efficiency $\epsilon_{\rm rad}$. There is a critical value (upper limit) $\dot{m}_{\rm crit}\sim\alpha^{2}\sim(0.01-0.1)$ for existing an ADAF (Narayan and Yi 1995). Given that the Bondi accretion rate is well below the Eddington rate (see Sect. 4.1), the BH accretion in M 87 is certainly in the RIAF regime, where the disk is optically thin and hot, and thus becomes geometrically thick (the ion is heated up to the virial temperature $T_{\rm i}\lesssim 10^{12}\,{\rm K}$, while the electron temperature is much cooler $T_{\rm e}\sim 10^{9-11}\,{\rm K}$: Nakamura et al 1997; Quataert and Gruzinov 1999).

In the quiescent state (QS) far below the critical mass accretion rate (say $\dot{m}<10^{-3}$), an ADAF-type accretion flow occupies all the way down to the innermost stable circular orbit (ISCO) around the BH. However, there is a state transition from an ADAF to a radiatively efficient (luminous), geometrically thin, standard accretion disk (Shakura and Sunyaev 1973; Novikov and Thorne 1973). A large transition radius $R_{\rm tr}\gtrsim 10^{4}\,r_{\rm S}$ is expected in the QS such as M 87 and Sgr A^∗ (Narayan and McClintock 2008). In the low state (LS) where $\dot{m}$ is higher than that in QS (i.e., $\dot{m}\sim 10^{-3}-10^{-1.5}$), the radiative efficiency becomes larger, and the transition radius gets smaller. In the intermediate state (IS), where $\dot{m}\sim(10^{-1.5}-10^{-1})\lesssim\dot{m}_{\rm crit}$, the radiative efficiency is as high as $\sim 0.1$ and the transition radius approaches the ISCO. Finally, the state transition is completed in the high state (HS), where $\dot{m}\gtrsim\dot{m}_{\rm crit}$, and the accretion flow changes completely to a thin disk down to the ISCO. This system is in the thermal state so that its spectrum is well-described by a multi-color disk model (see also Esin et al 1997; Narayan and McClintock 2008; Ho 2009, for more details).

Among the jetted RLAGNs, low power radio galaxies such as FR-I sources (e.g., M 87) can fall into the QS. High power radio galaxies such as FR-II sources (e.g., Cygnus A) can be in the LS. Ghisellini et al (2010) studied all blazars of known redshift detected by the Fermi satellite and provided a possible ‘split’ between BLOs and FSRQs in the accretion mode, which occurs at $\dot{m}\sim\dot{m}_{\rm crit}$ (see also Xu et al 2009; Chen et al 2015). Based on these observational facts, we can categorize BLOs in the LS and the IS, while FSRQs in the HS. Furthermore, some Fermi FSRQs may be in the very high state (VHS) where $\dot{m}\gtrsim 1$. This is the so-called the ‘super-Eddington’ accretion regime in the optically thick RIAF or slim disk (Begelman and Meier 1982; Abramowicz et al 1988). Also, NLSy1s are considered to be in the VHS (Collin and Kawaguchi 2004); transient and highly variable NLSy1s may indicate that such high-amplitude variability is linked with a transition between a standard disk and a slim disk as a result of the thermal instability, with $3\lesssim\dot{m}\lesssim 300$ (Kawaguchi 2003).

The magnetic field that an isolated BH can possess is weak: if there is a strong magnetic field, the mass accretion is prevented by its magnetic pressure and torque. Blandford and Znajek (1977) estimated the critical field strength for a maximally spinning BH as $\approx 2.5\times 10^{5}\ {\rm G}\ (M_{\rm BH}/M_{\odot})^{-1/2}$. On the other hand, the field strength brought in by accreting matter can be much stronger. Thus, an intrinsic magnetic field of the BH can be considered as astrophysically unimportant.

GRMHD simulations reveal that the poloidal magnetic flux in the accretion flow is advected inward until the ram pressure of accreting material ($\lesssim\rho c^{2}=GM_{\rm BH}\rho/r_{\rm g}$, shown as $F_{\rm g}$ in Fig. 3) is balanced with the magnetic pressure ($B^{2}/8\pi$, where $B$ is the strength of magnetic field lines, shown as $F_{\rm m}$ in Fig. 3). The accumulated poloidal magnetic flux onto the event horizon threading (one hemisphere of) the BH ($\Phi_{\rm BH}$) is a key parameter that essentially characterizes the accretion flow properties. This is given in the following form (Tchekhovskoy et al 2011; Yuan and Narayan 2014):

When the dimensionless parameter $\phi_{\rm BH}\gtrsim 50$ (in Gaussian units), a strong magnetic field pushes aside the accretion flow and tends to produce outflows. This is known as the Magnetically Arrested Disk (MAD) state (Igumenshchev et al 2003; Narayan et al 2003; Tchekhovskoy et al 2011). On the other hand, the case $\phi_{\rm BH}\sim a~{}few$ corresponds to the Standard Accretion and Normal Evolution (SANE) state, where weak and turbulent magnetic fields dominate the accretion flow (Narayan et al 2012).

The jet launching mechanism in AGNs is a longstanding challenge in astrophysics. Electromagnetic radiation emitted from AGN jets spanning from radio to very-high-energy (VHE) $\gamma$ rays have improved our understanding of AGN jets. The total power of the jet ($L_{\rm j}$) is the sum of the several components as described by

where $L_{\rm rad}$, $L_{\rm p}$, $L_{\rm e^{\pm}}$, and $L_{\rm Poy}$, are the radiation power, the proton power, the electron-positron pair power, and the Poynting power of the jet, respectively. (e.g., Sikora et al 2020).

If the Blandford and Znajek (BZ) process (Blandford and Znajek 1977), where the electromagnetic energy extraction from a rotating BH is operated as a jet launching mechanism, $L_{\rm Poy}=L_{\rm BZ}$ holds at the jet base and BZ power ($L_{\rm BZ}$) can be given by

where $f(\Omega_{\rm H})=1+1.38(\Omega_{\rm H}r_{\rm g}/c)^{2}-9.2(\Omega_{\rm H}r_{\rm g}/c)^{4}$, $\kappa$ is the geometrical factor, and $\Omega_{F}=\Omega_{H}/2$ is assumed (see Tchekhovskoy et al 2010, 2011, for details).

Observed AGN luminosities form two separated sequences on the radio-loudness–Eddington-ratio plane, suggesting larger BH spins for radio-loud sources (Sikora et al 2007). Zamaninasab et al (2014) suggest dynamically important magnetic fields near accreting BH based on the estimation of $\Phi_{\rm BH}$. They seem to support the idea of the BZ process at work and thus, most of theoretical studies assume that magnetic fields play a key role in jet formation. However, it is not obvious whether observed electromagnetic signals indeed suggest that the magnetic field plays a major role in the dynamics or not. This assertion should be verified through comprehensive observations, and M 87 can provide us an excellent platform to explore this (see Sect. 5.10).

The $L_{\rm rad}$ is directly obtained from observational data and they are mostly dominated by the power of synchrotron emission ($L_{\rm syn}$) and inverse-Compton emission ($L_{\rm IC}$) (Kataoka and Stawarz 2005). Based on the works of (Sikora et al 2007; Ghisellini et al 2014) the observed bolometric luminosity is known to be in a wide range of $10^{40}~{}{\rm erg~{}s^{-1}}\lesssim L_{\rm rad}\lesssim 10^{48}~{}{\rm erg~{}s^{-1}}$. On the other hand, constraining the plasma composition and estimating $L_{\rm p}$, and $L_{\rm e^{\pm}}$ in AGN jets remains challenging and is still a subject of debate since the observed emissions predominantly originated from only nonthermal $e^{\pm}$ pair (Wardle et al 1998; Sikora and Madejski 2000; Kataoka et al 2008; Kino et al 2012; Kawakatu et al 2016; Sikora et al 2020).

Although the estimation of $L_{\rm j}$ contains a certain degree of ambiguity due to the reason mentioned above, the jet power of FSRQs can maximally reach up to $L_{\rm j}\sim 10^{49}~{}{\rm erg~{}s^{-1}}$, while that of radio galaxies can reach as low as $L_{\rm j}\sim 10^{42}~{}{\rm erg~{}s^{-1}}$ (see for example Rawlings and Saunders 1991; Daly 2019; Chen et al 2023b, and references therein).

The mass of the central SMBH ($M_{\rm BH}$) of M 87 has been estimated using three different approaches: stellar dynamics modeling (e.g., Gebhardt et al 2011), gas dynamical modeling (e.g., Harms et al 1994; Macchetto et al 1997; Walsh et al 2013), and BH shadow measurements (Event Horizon Telescope Collaboration et al 2019a, f; Event Horizon Telescope Collaboration et al 2024).

Stellar dynamical measurements involve assessing the central rise in stellar velocity dispersion towards the center of the galaxy. This method began with the assumption of a simple Gaussian distribution of stellar dynamics, then it has evolved to estimate the mass-to-distance ratio, $M/D$, with high accuracy through the development of sophisticated modeling techniques utilizing the full line-of-sight stellar velocity distributions achieved with high angular resolution spectroscopic observations. As a byproduct, it provides the mass-to-light ratio, $M/L$, of the stellar populations at its center, which can be used to validate measurements by comparing them with independent estimates derived from population synthesis. Utilizing these approaches, $M_{\rm BH}$ of M 87 was estimated to be $M_{\rm BH}=(6.6\pm 0.4)\times 10^{9}\,M_{\odot}$ for an assumed distance of 17.9 Mpc (corresponding to $6.2\times 10^{9}\,M_{\odot}$ for 16.8 Mpc) (Gebhardt et al 2011). Recently, new estimates on $M_{\rm BH}$ have been proposed by adopting new triaxial orbit models, resulting in a relatively smaller mass of ($5.3_{-0.22}^{+0.37})\times 10^{9}\,M_{\odot}$ (Liepold et al 2023). On the other hand, a larger mass of $8.7\times 10^{9}\,M_{\odot}$ has been suggested by new measurements using the Canada–France–Hawaii Telescope (CFHT) and the Very Large Telescope (VLT) adaptive optics observations (Simon et al 2024). However they also estimated a smaller mass of $5.5\times 10^{9}\,M_{\odot}$ depending on the adopted stellar density profile (Simon et al 2024).

Gas dynamical estimates rely on modeling the dynamics of the H$\alpha$ and [N II] emission lines from the central gas disk. In early 1990s, Harms et al (1994) detected a spatially resolved ionized gas disk with the HST. They identified an asymmetry in the velocity field and estimated $M_{\rm BH}=(2.6\pm 0.8)\times 10^{9}\,M_{\odot}$. Following this, Macchetto et al (1997) revealed the rotation curve of the ionized gas disk surrounding the nucleus. The velocities measured from the emission lines were consistent with the hypothesis that the ionized gas disk undergoes Keplerian rotation, leading to $M_{\rm BH}=(3.6\pm 1.0)\times 10^{9}\,M_{\odot}$. Later, Walsh et al (2013) further refined the analysis and derived $M_{\rm BH}=(3.3\pm 0.8)\times 10^{9}\,M_{\odot}$ based on comprehensive gas-dynamical modeling with new HST data acquired with the Space Telescope Imaging Spectrograph (STIS). Historically, these $M_{\rm BH}$ estimates using gas dynamics have shown excellent agreement with each other, whereas there has been a factor of $\sim$2 discrepancy with the $M_{\rm BH}$ derived from the stellar dynamics.

While the above two methods are based on the orbital motion of gas or stars relatively far from the central BH, it is also possible to estimate $M_{\rm BH}$ by directly measuring the size of the BH shadow using ultra-high-angular-resolution VLBI observations at millimeter wavelengths. According to GR, the apparent size of the BH shadow is expected to be approximately 5 times the Schwarzschild radius (e.g., Event Horizon Telescope Collaboration et al 2019f). The first observational hint for such a gravitationally-lensed shadow feature was obtained for Sgr A* and M 87, with proto-EHT experiments at 1.3 mm using three stations in the United States (Doeleman et al 2008, 2012). The first-ever BH shadow image of M 87 obtained with EHT discovered an asymmetric ring-like structure with a diameter of 42 $\pm$ 3 $\mu$as (Event Horizon Telescope Collaboration et al 2019a, d, Fig. 4 left).

Comparing the observed M 87 EHT images with a large GRMHD simulation library, they concluded that the observed ring-like structure is in good agreement with the shadow of a spinning Kerr BH (Event Horizon Telescope Collaboration et al 2019e). By calibrating the observed ring diameter with the predicted ring diameter based on the simulations, they derived the gravitational radius and estimated $M_{\rm BH}$ to be $M_{\rm BH}=(6.5\pm 0.7)\times 10^{9}\,M_{\odot}$ assuming a distance of 16.8 Mpc (Event Horizon Telescope Collaboration et al 2019f). The first EHT results based on the 2017 observations have been confirmed with another independent EHT dataset taken in 2018 (Event Horizon Telescope Collaboration et al 2024, Fig. 4 right), providing strong supporting evidence that the ring seen in the EHT images represents the BH shadow predicted by GR.

The $M_{\rm BH}$ estimated from the EHT data is in good agreement with that derived from the stellar dynamics, while it deviates from the ones estimated from the gas dynamics with a confidence level of 99 $\%$ (Event Horizon Telescope Collaboration et al 2019f). This discrepancy may be attributed to the inclusion of internal gas motions within the disk, along with adopting slight changes in the inclination angle of the disk (Jeter et al 2019; Jeter and Broderick 2021). Therefore a further study is warranted, particularly using future high-resolution optical instruments.

The spin of a BH is a fundamental parameter and crucial for describing its nature. In astrophysical contexts, the BH spin emerges as a compelling explanation for the extraordinary power exhibited by jets. Take, for example, the case of M 87, where the estimated mass accretion rate onto the SMBH ranges from 10^-3 to 10^-4 $M_{\odot}\,{\rm yr}^{-1}$. This translates to a mass accretion power $\dot{M}c^{2}$ of $\sim$10^43-44 erg s^-1 (as discussed in Sect. 4.3). However, the expected kinematic power of the jets falls in the range of $\sim$10^43-45 erg s^-1 (as outlined in Sect. 5.11). Hence, it requires an energy conversion efficiency of $\gtrsim 100\%$, implying that the simple mass accretion power alone may be insufficient to explain the jet power. It is postulated that the subtraction of the rotating energy of the BH can compensate for this difference.

The direct imaging of SMBH shadow provides another important clue regarding the nature of BH beyond its mass. The observed brightness asymmetry of the 1.3 mm ring is likely caused by the effect of Doppler beaming associated with the frame dragging effect of a spinning Kerr BH. Since the observed ring appears brighter in the southern part of the emission, it suggests that the BH spin vector is oriented away from Earth if the BH spin axis aligns with the large-scale jet (Event Horizon Telescope Collaboration et al 2019f). The bright spot observed in 2018 shifted by 30 degrees relative to that observed in 2017, and the asymmetry observed in 2018 is more consistent with the orientation of the large-scale jet (Event Horizon Telescope Collaboration et al 2024). These observed features indicate the presence of a spinning BH at the center of M 87, although the value of the spin has not yet been tightly constrained in these observations. Further strong evidence of the spinning black hole will be provided by precise measurement of the photon ring with next generation VLBI observations (see Sect. 6.2), and the next step is to determine the direction of the spin relative to the rotation of the accreting matter (prograde/retrograde).

As described in Sect. 2.3, the nature of accretion flows onto SMBH can be probed by observing the X-ray hot gas atmosphere. The density and temperature profiles of the X-ray hot gas in M 87 on kpc scales were extensively investigated using Chandra X-ray data (Di Matteo et al 2003; Russell et al 2018). Di Matteo et al (2003) inferred the decrease in the gas temperature from 10 kpc down to the central 1 kpc. Within this distance range, they observed the density profile to flatten inside the central 2 kpc, while it follows $r^{-1}$ beyond 2 kpc up to 10 kpc (see Fig. 5). With these observed profiles, they calculated the Bondi radius to be $\sim$150 pc for a BH mass of $3\times 10^{9}$ $M_{\odot}$ and estimated the Bondi accretion rate to be $\dot{M}_{\rm B}\sim 0.1\,M_{\odot}\,\rm{yr}^{-1}$. The temperature and density profiles were further refined by Russell et al (2018) through a stacking analysis of twelve 5-ks Chandra data. By adopting the refined temperature and density profile together with a newly estimated $M_{\rm BH}$ of 6.5 $\times$ 10⁹ $M_{\odot}$ (Event Horizon Telescope Collaboration et al 2019a, f), we estimate the Bondi radius and Bondi accretion rate to be

and

Hence, the comparison with $\dot{M}_{\rm Edd}$ introduced in Eq.(5) ($\simeq 140\,M_{\odot}\,{\rm yr^{-1}}$ for $\epsilon_{\rm rad}=0.1$ and $M_{\rm BH}=6.5\times 10^{9}\,M_{\odot}$) indicates that the accretion onto the M 87 nucleus is significantly sub-Eddington, as also inferred from the underluminous core X-ray luminosity ($\sim 7\times 10^{40}$ erg s^-1; Di Matteo et al 2003). It is notable that Russell et al (2018) found no evidence of gas temperature increase within a radius of 0.25 kpc, as expected in the classical Bondi-type accretion flows, while a further increase in density profile was observed within 0.3 kpc. This led them to suggest that the classical Bondi accretion flow may not be accurately achieved, and that a substantial decrease in the mass accretion rate onto the central SMBH within $r_{\rm B}$.

In addition to probing the hot plasma through X-ray observations, efforts have been made to unveil the cold accretion flow at similar or slightly smaller scales using radio interferometry data. Tan et al (2008) utilized the Submillimeter Array (SMA) at 230 GHz to detect CO (J = 2–1) line emission, aiming to identify the presence of cold molecular gas in a thin disk around the M 87 BH. Although they detected 230 GHz continuum emission from the nucleus and several knots at arcsecond scals, they estimated a conservative upper limit on the mass of molecular gas to be $\sim$ 8 $\times$ 10⁶ $M_{\odot}$ within 100 pc of the central BH. Subsequently, Simionescu et al (2018) detected extended CO (2–1) line emission approximately 40 arcseconds away, corresponding to 3.4 kpc in the southeast direction from the nucleus, using the Atacama Large Millimeter/submillimeter Array (ALMA). They derived the corresponding molecular gas mass to be (4.7 $\pm$ 0.4) $\times$ 10⁵ $M_{\odot}$. Following this, Li (2022) conducted a systematic survey utilizing ALMA archival data. However, they did not find conclusive evidence of CO line emission in the vicinity of the M 87 nucleus.

It is difficult for the current X-ray instruments to spatially resolve the scales well below $r_{\rm B}$. To explore the nature of accretion flows within $r_{\rm B}$ scales, a powerful approach would be the measurement of polarization and Faraday rotation measure (RM) at mm/submm wavelengths since the bulk of mm/submm synchrotron emission from ADAF/RIAF-type accretion flows is likely originated in the close vicinity of the central SMBH (e.g., Yuan and Narayan 2014). This approach was originally developed for the nearest SMBH Sgr A*, which exhibits highly linear polarization with significant RM (e.g., Aitken et al 2000; Bower et al 2003; Marrone et al 2006, 2007). By attributing the observed RM to the magnetized plasma associated with the accretion flow, the accretion rate of Sgr A^∗ at small radii was constrained to be $\dot{M}\sim 10^{-8}\,M_{\odot}\,\rm{yr}^{-1}$ depending on the adopted accretion models (Marrone et al 2006, 2007).

Following the success in constraining the $\dot{M}$ of Sgr A*, Kuo et al (2014) applied a similar method to M 87 using 1.3 mm polarimetric data taken by the SMA. They obtained an upper limit of $\lvert{\rm RM}\rvert$ (the magnitude of RM) to be $7.5\times 10^{5}\,{\rm rad\,m}^{-2}$. Adopting a simple analytic model of the accretion flow similar to the case of Sgr A*, the $\dot{M}$ onto the M 87 SMBH was constrained to be $<9\times 10^{-4}\,M_{\odot}\,\rm{yr}^{-1}$. This indicates a substantial decrease in $\dot{M}$ compared to the Bondi accretion rate. This is in good agreement with expectations from the RIAF models, ruling out the classcial ADAF.

Large RM values towards the M 87 nucleus were later confirmed by Goddi et al (2021) using multi-epoch ALMA data at 3 and 1.3 mm. Thanks to the improved sensitivity of ALMA, they further detected significant time variation of RM both in magnitude and sign. This suggests a more complicated origin of the Faraday rotation rather than a simple accretion flow, making a more accurate estimate of $\dot{M}$ rather challenging.

To look into the structure inside $r_{\rm B}$ more directly, high-resolution VLBI observations are required. Using the Very Long Baseline Array (VLBA) at cm wavelengths, there were several attempts to explore the pc-scale polarimetric and RM properties of M 87 (e.g., Junor et al 2001; Zavala and Taylor 2002; Park et al 2019b, see also Sect. 5.5). Most notably, Park et al (2019b) revealed a spatially-resolved RM profile inside $r_{\rm B}$ and found that $\lvert{\rm RM}\rvert$ decreases with distance from 10 mas to 350 mas (corresponding to deprojected distances 5000–200000 $r_{\rm S}$) away from the nucleus. They indicated that the observed radial slope of RM at these scales is reproduced by a gas density profile of $\rho(r)\propto r^{-1}$ for the Faraday screen. Such a density profile is expected for RIAF-type hot accretion flows with substantial winds or mass loss (ADiabatic Inflow-Outflow Solutions (ADIOS); Blandford and Begelman 1999), which is consistent with the significant reduction of $\dot{M}$ near SMBH inferred from the mm/submm interferometric RM studies.

We note that optical observations also provide additional constraints. Besides constraining the SMBH mass, the HST data that detected the nuclear gas disk can be used to further infer the accretion rate associated with the gas disk. As described in Sect. 3.1, the mass estimated from the gas dynamics differs from those from the stellar dynamics and the EHT. The discrepancy could stem from inappropriate assumptions about the dynamics of the gas disk, as discussed in detail by Kormendy and Ho (2013). Jeter et al (2019) refined the gas-dynamical model by adopting a sub-Keplerian disk velocity together with a different viewing angle and demonstrated that the refined model can explain the larger BH mass. By using the refined model with conservative estimates of the cold gas at 100 pc scale by Tan et al (2008), they estimated the mass accretion rate associated with the gas disk to be $\sim$$2\times 10^{-3}\,M_{\odot}\,{\rm yr}^{-1}$.

The accretion flow structure near event-horizon scales has become accessible through direct imaging with mm-VLBI. As described in Sect. 3, the EHT observations at 1.3 mm unveiled a ring-like structure in the core of M 87 (Event Horizon Telescope Collaboration et al 2019a, d, f). By comparing the EHT images with an extensive library of GRMHD and GR radiative transfer (GRRT) simulated images, the observed ring-like structure was found to be in good agreement with the emission associated with the accretion flow (Event Horizon Telescope Collaboration et al 2019e). Further constraints on the properties of the 1.3 mm ring-like structure were obtained based on the subsequent analysis of linear-polarization data (Event Horizon Telescope Collaboration et al 2021a, b). It was reported that the observed electric vector polarization angle (EVPA) distributions predominantly describe a nearly azimuthal pattern, suggesting that the poloidal magnetic field component is more dominant than the toroidal one at these scales. An extensive comparison of the EHT polarization images with a large GRMHD simulation library indicated that MAD is the favored state of the accretion flow, rather than SANE ones, suggesting the presence of dynamically important magnetic fields near SMBH. The recent further detection of circular-polarization signals in the EHT 2017 data additionally reinforces the MAD scenario (Event Horizon Telescope Collaboration 2023). Based on the GRMHD models that are in good agreement with the observations, the mass accretion rate at the EHT scales was estimated to be $\dot{M}\sim(3-20)\times 10^{-4}\,M_{\odot}\,{\rm yr}^{-1}$ (Event Horizon Telescope Collaboration et al 2021b). Analysis of the EHT total-intensity and polarization images also provided horizon-scale estimates on other physical quantities such as electron density ($n_{\rm e}\sim 10^{4}\text{--}10^{7}\,{\rm cm}^{-3}$), electron temperature ($T_{\rm e}\sim(1\text{--}12)\times 10^{10}\,\text{K}$), and magnetic field strength ($B\sim(1-30)$ G).

Following the EHT imaging of the BH shadow at 1.3 mm, recent Global Millimeter VLBI Array (GMVA) observations connected to ALMA and the Greenland Telescope (GLT) have revealed another ring-like image in the M 87 nucleus at a wavelength of 3.5 mm (Lu et al 2023, see also Fig. 6). The diameter of the 3.5 mm ring is approximately 8.4 times the Schwarzschild radius, which is 50% larger than the diameter of the 1.3 mm ring. Additionally, the inner edge diameter of the 3.5 mm ring is also larger than the diameter of the EHT ring. The larger and thicker ring size at 3.5 mm is likely caused by the opacity effect of synchrotron emission at a lower frequency. A comparison of the GMVA image with GRMHD simulations in Lu et al (2023) provides a strong indication that the observed ring-like structure at 3.5 mm is dominated by the accretion flow.

As a summary of this section, in Fig. 7 we show a plot compiling various estimates of $\dot{M}$ in the literature, measured at various scales. Although the statistics are still very limited, a collection of these estimates reveals a substantial decrease in mass supply towards the innermost regions with respect to the Bondi accretion rate. This is clear evidence that the classical ADAF is no longer valid as the model for the M 87 accretion flows, while the suggested radial evolution of $\dot{M}$ is consistent with with ADIOS-like solutions.

While the EHT 230 GHz observations resolved the shadow of SMBH at the jet base, it is worth noting the nature of the ‘radio core’ seen at the lower radio frequencies. The radio core, a bright compact feature at the apparent origin of an AGN jet seen in a VLBI image, is most likely either a synchrotron-self-absorbed opaque surface of the jet base at a given frequency (e.g., Konigl 1981; Lobanov 1998) or a standing shock feature (e.g., Daly and Marscher 1988; Marscher et al 2008). In the case of M 87, multi-frequency VLBA astrometric observations measured a core-shift between 2 and 88 GHz with a frequency dependence of $\sim\nu^{-(0.9-1.0)}$ (Hada et al 2011; Jiang et al 2021), in agreement with the M 87 radio core at these frequencies being a photosphere of the jet base within a few tens $r_{\rm S}$ from the central BH. On the other hand, as introduced in Section  4, the recent GMVA+ALMA+GLT 86 GHz observations have spatially resolved the radio core into a ring-like structure, with emission likely dominated by the inner part of accretion flows (Lu et al 2023). Thus at $\sim$3 mm, the emission from both the accretion flows and the jet base could be blended when the radio core is imaged with a lower angular resolution.

At cm wavelengths. a handful of ultra-high-resolution space-VLBI observations resolved the radio core region into highly complex morphology (Dodson et al 2006; Asada et al 2016; Kim et al 2023). Especially at lower frequencies (1.6/5 GHz), space-VLBI images revealed diffuse substructures surrounding the compact core perpendicular to the jet axis (Dodson et al 2006; Asada et al 2016). The nature of this low-level emission has not yet been explored in detail, but it may hint at the presence of transverse stratification in the base of jet/outflows.

The jet of M 87 is predominantly one-sided on both kpc and pc scales, indicating its strong Doppler-boosted nature. This suggests that the viewing angle of the jet is relatively small or/and the jet speed is highly relativistic at these scales. On parsec scales, hints of the counter jet were observed in some early VLBI images (e.g., Reid et al 1989; Junor and Biretta 1995), but it did not receive much attention until Ly et al (2004) claimed the first detection of a faint counter jet at the eastern side of the radio core, although the possibility of a calibration artifact was not entirely ruled out.

The (sub)pc-scale counter jet of M 87 is now routinely observed in high-quality VLBI images at various frequencies (Fig.10), and careful imaging tests by multiple teams demonstrate that this feature is very likely real (e.g., Kovalev et al 2007; Ly et al 2007; Kim et al 2018b; Walker et al 2018; Nikonov et al 2023). Interestingly, Kovalev et al (2007) and Walker et al (2018) reported that the counter jet was also edge-brightened similar to the main jet, and a similar structure was also identified in more recent study by Nikonov et al (2023) (Fig. 10). The counter jet appears to extend $\gtrsim$3 mas from the core at 15 GHz, but looks even more extend if we take a look at the high-sensitivity VLBI image at 8 GHz presented in Nikonov et al (2023). The kinematics of the counter jet are much less constrained than that of the main jet, but independent measurements by different groups (Kovalev et al 2007; Ly et al 2007; Hada et al 2016; Mertens et al 2016) consistently reported substantially slow apparent speeds ($\sim$0.01–0.17 $c$).

The jet-to-counter-jet brightness ratio depends on various factors such as observing frequency, core distance, possible variability, and whether the ratio is calculated using integrated fluxes or peak intensity. However, a typical range would be of the order of 10 at $\sim$1 mas from the core. Combining this with proper motion measurement of the jet and counter jet, a conservative estimate of the jet viewing angle results in $\theta_{\rm view}$$\sim$13–40^∘. A tighter constraint of $\theta_{\rm view}$$\sim$14–20^∘ is obtained by considering more detailed jet velocity field (Mertens et al 2016), consistent with the angle suggested from the HST-1 superluminal motion (Sect. 5.7).

Since the first identification of polarized emission from M 87 by Baade (1956), the polarimetric properties of this jet have been extensively investigated across various scales and wavelengths. Polarimetry serves as a powerful tool for probing the structure of magnetic field, the internal jet structure as well as interactions between the jet and its surrounding environment.

On kpc scales, the polarization of the M 87 jet has been studied in great detail since 1980s using VLA in cm radio bands (e.g., Owen et al 1989, 1990) and later also using HST in optical bands and ALMA in mm bands (e.g., Capetti et al 1997; Perlman et al 1999; Avachat et al 2016; Goddi et al 2021). Polarized emission is detected throughout jet regions on kpc scales, with enhanced signals predominantly associated with the bright knot regions. While the observed magnetic field position angles (MFPAs; formally $90^{\circ}$ rotated from the measured EVPAs) on kpc scales are largely parallel to the jet axis in both radio/optical bands, only the optical MFPA vectors tend to become perpendicular to the jet at the upsteam end of some bright knots, suggesting a multi-layered structure in the jet (Perlman et al 1999). More recently, analyzing a wideband (4–18 GHz) Jansky VLA full-polarimetric data, Pasetto et al (2021) revealed extremely detailed radio polarization and Faraday RM distributions that support the presence of a helical magnetic field on kpc scales. See also Sect. 5.8 for a review on the kpc-scale jet.

In contrast to the rich polarization structures seen at kpc and horizon scales (see Sect. 4.3), still much less is known about the polarimetric properties on intermediate ($\sim$0.1–100 pc) scales, as the pc-scale jet of M 87 is largely unpolarized at radio wavelengths. Pioneering studies were conducted by Junor et al (2001) and Zavala and Taylor (2002) using VLBA at cm wavelengths (5–15 GHz), where they detected patchy polarization features in the jet portion around $\sim$20 mas from the core with considerably large RM magnitudes (thousands to $\sim$$10^{4}$ rad m^-2). Weak patchy polarization features were also detected closer to the core at 86 GHz (Hada et al 2016). A significant advance was made by Park et al (2019b), where they revealed more detailed RM distributions along the pc-scale jet by revisiting multi-frequency polarimetric VLBA data at 2–8 GHz. They found that the RM magnitude gradually decreases with distance from one parsec scales to just before HST-1. The observed slope of RM distributions was reproduced with the scenario that the Faraday screen is the winds from the hot accretion flows, which may be relevant to the observed collimation and acceleration of the M87 jet on the same scales. More recently, Nikonov et al (2023) have reported the presence of an RM gradient also in the direction perpendicular to the jet based on their high-sensitivity VLBI images at 8 and 15 GHz.

Between the horizon scales and parsec scales, the polarization structure of the radio core at subpc scales remains controversial. Walker et al (2018) and Kravchenko et al (2020) reported complex EVPA patterns that appeared to be wrapped around the 43 GHz core. However, a later analysis by Park et al (2021a) with an improved calibration method did not find such a complex feature but detected a simple compact polarization component with its peak coinciding with the total intensity peak. To obtain a definitive answer and bridge the magnetic field structure between horizon and parsec scales, higher sensitivity VLBI observations at 43/86 GHz combined with a refined polarization calibration technique would be required.

As mentioned in Sect. 5.1, the pc-scale jet of M 87 is widely known to show a limb-brightened structure with a parabolic shape. However, when imaged at a higher resolution or/and sensitivity, another streamline appears to emerge in the middle of the jet (Fig. 11). Such a ‘spine’ component was clearly detected in high-quality VLBI images at cm wavelengths (Asada et al 2016; Hada 2017; Tazaki et al 2023; Nikonov et al 2023) and also in stacked VLBI images at 43/86 GHz (Walker et al 2018; Kim et al 2018b). The central spine emission appears to be substantially narrower than the whole jet width. Recently, a new GMVA 3 mm M87 image connected to ALMA and GLT has clearly detected a triple-ridge jet structure even closer to the core (Lu et al 2023). These images obtained at various scales/frequencies indicate that such a ‘spine-sheath’ structure of this jet is maintained over a wide range of distances from the BH vicinity to further out.

The nature of the observed triple-ridge profile remains a matter of debate, with some potential scenarios proposed to explain the origin. One possibility is that the central ridge is actually part of the outer sheath associated with the same layer as the northern/southern limbs, while the true spine emission near the jet axis is less dominant due to either a low emissivity or deboosting (Mertens et al 2016; Walker et al 2018). solely reproduced by a BH-powered jet (Asada et al 2016; Orihara2019), with different streamlines or beaming factors. Another possibility is that the central ridge is associated with the true spine within the jet interior, which may be anchored in the central rotating BH (Asada et al 2016; Sob’yanin 2017; Ogihara et al 2019).

At a deprojected distance of $\sim$200 pc from the nucleus, the M 87 jet harbors a notable feature known as ‘HST-1’. This feature was initially discovered with HST by Biretta et al (1999), and has been attracting special attention from the high-energy astrophysical community for several reasons.

First, HST-1 contains highly superluminal components of $\sim$(4–6) $c$ that are observed across radio, optical and X-ray bands (e.g., Biretta et al 1999; Cheung et al 2007; Giroletti et al 2012; Snios et al 2019). The observed apparent speeds appear to be the highest in the entire velocity field of the M87 jet (Meyer et al 2013; Asada et al 2014, see also Fig. 8 bottom). This suggests that HST-1 represents a maximally Doppler-boosted region in the M 87 jet, and also tightly constrain the jet viewing angle to be $\theta_{\rm view}\lesssim 19^{\circ}$. Second, as elaborated further in Sect. 5.9, HST-1 experienced a notable MWL outburst around 2005 (Harris et al 2006), followed by the ejection of superluminal components (Cheung et al 2007, Fig. 12 top). Third, HST-1 is located in the transition region of the jet geometry, where the jet shape changes from parabolic to conical. Additionally, there is a contraction of the jet cross section at HST-1 itself (Asada and Nakamura 2012), accompanied by highly enhanced polarization and RM (Park et al 2019b). The combination of these observational characteristics has led to an intriguing hypothesis that HST-1 represents a recollimation shock region at the end of ACZ in the M 87 jet, releasing significant energy in the form of high-energy flares. HST-1 could therefore be a possible counterpart of the unresolved core of distant blazars.

When viewed at the highest resolution, HST-1 exhibits intricate internal structure and time evolution (e.g., Chang et al 2010). It consists of multiple substructures with different velocities and trajectories. Long-term VLBI monitoring observations of HST-1 detected repeated ejection of new superluminal components from the upstream edge of the HST-1 complex (‘HST-1d’; Cheung et al 2007; Giroletti et al 2012; Hada et al 2014a, see also an updated plot in Fig. 12), suggesting that the upstream edge is where particles passing through are reenergized. After the ejection, while all the components are ejected from a similar location, their subsequent trajectories are quite different from one another, with most components following curved paths. This implies that the actual trajectories of these components are three-dimensional. The presence of helical motion in HST-1 is suggested by Chen et al (2011), where they detected a progressive rotation of HST-1’s EVPA using VLA data around the 2005 event. Unfortunately, the brightness of HST-1 has been continuously decreasing since the 2005 event across all wavelengths, rendering it very faint and challenging to observe as of 2024. Reenergization of HST-1 is eagarly awaited.

In this section, we review the estimations of some key quantities relevant to the strength of magnetic field based on actual observational data of M 87.

As reviewed in Sect. 2.5, it is difficult to put a limit on the total amount of nonthermal electrons (and positrons) because information can only be obtained directly from the jet’s electromagnetic observations for a portion of nonthermal electrons (and positrons). Therefore, there have been no significant recent progress in the estimation of $L_{\rm j}$ for M87, and we will rely on the estimates made in relatively past studies. In the following, we proceed with the review in order of increasing spatial scale.

On several 10 kpc scales, $L_{\rm j}$ is estimated from the energy required to inflate the observed X-ray cavity and radio lobe structures. The previous works suggest $L_{\rm j}\sim 10^{43-44}~{}{\rm erg~{}s^{-1}}$ (Owen et al 2000; Young et al 2002; Allen et al 2006; Rafferty et al 2006). The spherical shock model provides an estimate of the mean power of the shock outburst as $2\times 10^{43}~{}{\rm erg~{}s^{-1}}$ averaged over the past $10^{7}$ year on $\sim 50$ kpc scale (Forman et al 2005). An application of synchrotron aging model against the Low-Frequency Array (LOFAR) observation of the radio lobe suggests an age of $\sim$40 Myr, which leads to $L_{\rm j}\sim(6-10)\times 10^{44}~{}{\rm erg~{}s^{-1}}$ (de Gasperin et al 2012).

At a distance of about 1 kpc from the central BH, Knot A is known to be a bright optical feature in the jet downstream side. Bicknell and Begelman (1996) made an estimation of $L_{\rm j}$ using the Knot A. Interpreting knot A as an oblique shock within the jet results in a jet power estimate of $(1-3)\times 10^{44}~{}{\rm erg~{}s^{-1}}$. The estimated $L_{\rm j}$ is consistent with the above-mentioned jet power required to feed the radio lobe. In addition, the search for a steady and extended $\gamma$-ray signal around M 87 can constrain the upper limit of cosmic-ray energy density, as discussed in a recent study by the H.E.S.S. collaboration, which addressed the upper limit of $L_{p}$ (H.E.S.S. Collaboration et al 2023).

By identifying HST-1 as a recollimation shock, Stawarz et al (2006) estimated $L_{\rm j}$ by the pressure balance between the ambient matter and HST-1. They derive a jet power of the order of $10^{44}~{}{\rm erg~{}s^{-1}}$. The estimated $L_{\rm j}$ is again consistent with the above-mentioned jet power required to feed the radio lobe.

In addition, it is worth noting that Broderick et al (2015) pointed out that the above-mentioned $L_{\rm j}$ implies a certain minimum $\dot{M}$ and if the central compact object in M 87 were not a BH but had a surface, the emission from the accreting matter onto the surface would exceed the observed fluxes in near-infrared and optical bands, suggesting the presence of an event horizon.

In Table 3 we summarize estimated values of various parameters and observational quantities for M 87, including those introduced in Sect. 2. The table reflects our current level of understanding about this source. One can see that some of the key quantities are already determined with sufficient accuracy, while some other parameters remain highly uncertain. It is also worth noting that, for some important physical quantities (e.g., the accretion rate, jet shape, jet velocity, magnetic field strength), we are starting to see their spatial evolution over many orders of magnitude in radial distance ($r$) or jet axial distance ($z$), owing to the accumulation of observational data at different scales and the progress of multi-scale studies over the past years. Therefore, we are now beginning to connect the physics between horizon scales and galactic scales, which is essential for a comprehensive understanding of this active galaxy as a whole system.

More than a quarter of a century has passed since the first dedicated workshop on the radio galaxy M 87 in Tegernsee, Germany in 1997. During the conference, some important questions were left unanswered (Blandford 1999): ‘Why is M 87 so dim?’; ’Does the Jet Power Come from the Disk or the Hole?’; ’How is the Jet Collimated and of What Is It Made?’; ’What’s Going on in Knot A?’. We believe that these fundamental quests have been extensively studied both observationally and theoretically since then. In addition, the second M 87 workshop⁴⁴4https://events.asiaa.sinica.edu.tw/workshop/20160523/ in Taipei, Taiwan, in 2016 drove momentum for the imaging of a BH shadow, along with preparations for the 100th anniversary of the discovery of astrophysical jets. Now, we know answers to the questions above and we know that there are SMBHs at the center of galaxies. This is, however, still a good time to be studying M 87, as Roger Blandford mentioned before.

Below we describe some of the selected remaining challenges or topics that should be explored in future studies of M 87.

• •

  BH spin: As summarized in Table 3, a growing number of studies suggest that the M 87 SMBH is likely spinning, whereas the exact value of the spin remains highly uncertain. Constraining the spin parameter is crucial for testing the spacetime properties near the SMBH and for testing the jet launching mechanism via BZ processes. To directly constrain the spin, one of the most promising ways would be to accurately measure the shape of a photon ring, where the ring shape is expected to deviate from a pure circular shape as the spin increases (up to $\pm$$\sim$5%) (Takahashi 2004; Johannsen and Psaltis 2010). The ring-like structure observed by the ground-based EHT baselines at 1.3 mm is likely still dominated by the emission from the gas surrounding the BH, while the emission from a lensed photon ring becomes significant only at baselines longer than $\sim$20 G$\lambda$ (Johnson et al 2020). To detect signals from such an uncontaminated photon ring, future higher-resolution mm/submm VLBI projects, such as the next-generation EHT (ngEHT⁵⁵5https://www.ngeht.org/), Black Hole Explorer (BHEX⁶⁶6https://www.blackholeexplorer.org/) with a space satellite and/or VLBI at even higher frequencies (e.g., GLT; Chen et al 2023a) will be promising.

• •

  BH-jet connection: While the EHT 2017/2018 data clearly detected the BH shadow, they failed to recover the surrounding jet emission mainly due to the limited baseline coverage of the array, leaving the exact connection between the shadow and jet launching still unresolved. Future EHT, ngEHT and/or BHEX observations at 230/345 GHz with significantly enhanced array capabilities hold promise for detecting both the shadow and the jet-base emission in the same image (e.g., Event Horizon Telescope Collaboration 2024), which is also crucial for constraining the BH spin as well as the jet-launching mechanism. This exploration will be further reinforced by complementing multi-frequency high-resolution images at 86 GHz and 43 GHz, through which one can trace a full evolution of initial jet formation from $\sim$1 $r_{\rm S}$ to $\sim$1000 $r_{\rm S}$.

• •

  Dynamical/Time-domain view: Thanks to the accumulation of multi-epoch, multi-wavelength, and multi-scale observational data, now we know that the structure of M 87 is variable at all spatial scales (from horizon to kpc) and time scales (from day to years). However, none of these observed variable features remain fully understood due to the complex nature of the underlying physics and the general lack of time-domain observational data. Time-domain/dynamical information provides much richer insights into the accretion and jet physics as well as the mechanisms of high-energy emission and particle acceleration processes. Continued multi-epoch, multi-instrument, multi-wavelength observations are highly encouraged to improve our knowledge of the time-domain properties of this object.

• •

  Mass transport processes: Dating back to the mid-1990s, ADAF (Narayan and Yi 1995) was proposed to account for SgrA* and the subsequent development of RIAF with substantial mass loss including the convection (Abramowicz et al 2002) and winds (Blandford and Begelman 1999) has advanced our understanding of the mass transport towards sub-Eddington BHs including M 87. A substantial decrease in $\dot{M}$ may be the norm, but we may not have a problem with this fact even if when considering jet production (Ghisellini et al 2005). When an accreting BH is in the MAD state (Narayan et al 2003), the spinning BH is able to power the jet (Blandford and Znajek 1977). However, the angular distribution of the accretion flow and winds remains unclear, so the net amount of mass transport onto the BH has not yet been confirmed at a quantitative level. The issue of mass and energy transport should be examined in multiple dimensions(Yuan and Narayan 2014; Yuan et al 2022). The fate of winds with significant mass loading is of particular interest, and further investigations with Faraday RM measurement are encouraged.

• •

  Jet structure: As reviewed in Sect. 5, a growing number of pc-scale VLBI images suggest that the transverse structure of the M 87 jet is likely more complicated than previously thought, as represented by the spine-sheath-like morphology. Additionally, these relativistic flows are likely surrounded by (invisible) non-relativistic outflows/winds. On kpc scales, high-quality VLA images confirm the double-helix morphology of the large-scale jet (a similar helical structure has recently been suggested also in the pc-scale jet; Nikonov et al 2023). Revealing the detailed internal structure of the M 87 jet is essential to better understand the formation, propagation and interaction with the external medium. To facilitate this study, further improvement in both angular resolution and image dynamic range for the extended jet is required. Future deep interferometric imaging observations with the next-generation VLA (ngVLA⁷⁷7https://ngvla.nrao.edu/) and the Square Kilometer Array (SKA⁸⁸8https://www.skao.int/) operated at cm wavelengths will enable us to conduct this kind of study.

• •

  Magnetic fields: Clarifying the geometry of magnetic field is crucial for understanding jet formation. On near-horizon scales, Chael et al (2023) pointed out that the sign of arg($\beta_{2}$) correlates with Poynting flux direction, and the EHT’s measurement of arg($\beta_{2}$)⁹⁹9The arg($\beta_{2}$) is the second azimuthal Fourier mode of linear polarization in a near-horizon images (Palumbo et al 2020). implies electromagnetic energy outflow. Chael et al (2023) found that (i) arg($\beta_{2}$) strongly depends on the BH spin, and (ii) arg($\beta_{2}$) shows a radial dependence, evolving rapidly near the inner shadow as fields are wound up. Thus whether this energy outflow is spin-powered or not remains inconclusive. To test whether the BZ mechanism really is at work or not in the M 87 jet, measuring arg($\beta_{2}$) closer to the event horizon is imperative. Higher-dynamic range or/and higher resolution polarization images with ngEHT/BHEX will be essential to clarify the radial dependence of arg($\beta_{2}$). On a parsec scale view of the M 87 jet, the linearly polarized emissions are weak and still detected only part of the jet (e.g., Zavala and Taylor 2002; Park et al 2019a). To explore the overall radial evolution of the magnetic field geometry and understand how it evolves from the BH to larger scales, future high-sensitivity VLBI and multi-wavelength polarimetric observations are highly anticipated.

• •

  Particle acceleration and high-energy emission: The origin of energetic particles emitting VHE $\gamma$-rays is of great interest, as described in Sect. 5.9. During the 2018 EHT campaign, flaring episodes constrained the size of the VHE $\gamma$-ray emitting region in SED modeling. However, a one-zone model cannot fully explain the SED, requiring an additional component. Two processes likely produce the energetic particles emitting VHE $\gamma$-rays: magnetic reconnection and shock acceleration. More sensitive VHE instruments such as the Cherenkov Telescope Array (Acharya et al 2013; Cherenkov Telescope Array Consortium et al 2019), may help identify the VHE $\gamma$-ray emitting region due to their better short-timescale sensitivity and spectral resolution. Investigating EVPA time evolution of the high-energy emission may also help distinguish these scenarios. Magnetic reconnection flares may cause drastic changes of EVPA, while shock-induced flares may not. Future X-ray polarization mission eXTP (Zhang et al 2019) could help distinguish these scenarios.

The extensive M 87 studies over the past decades have dramatically improved our understanding of various key properties of AGNs/SMBHs. Nevertheless, M 87 occupies only part of the whole AGN parameter space. To explore the diverse properties of AGNs and how M 87 fits into the context of a ‘unified picture’ of AGNs, comparison with other sources is essential.

First, the Galactic Center SgrA* is one of the most frequently compared objects with M 87 since they are the only two sources where the SMBH shadow is spatially resolved (Event Horizon Telescope Collaboration et al 2022a, b, c, d, e, f). Unlike M 87, SgrA* does not have a clear jet. This sharp contrast offers a unique opportunity to tackle the question, directly at the horizon scales, of why some SMBHs exhibit powerful jets while others do not. Interestingly, the EHT results of SgrA* favor the MAD state and high spin models, which might be expected to generate a jet efficiently like M 87. Other primary differences are that SgrA* has a much lower $M_{\rm BH}$ and $\dot{M}$ than M 87 both in absolute terms and when scaled by mass, which may also limit the jet power.

Second, while our understanding of weakly accreting SMBHs has greatly improved thanks to the intensive studies on M 87 (and SgrA*), the physics of accretion/ejection for other types of RLAGNs, especially those with highly accreting SMBHs, remains far less understood. Such sources may produce even more powerful jets, but they are typically much more distant than M 87. Recent ultra-high-resolution mm-VLBI and space-VLBI observations have started to resolve the inner jet regions within 100-1000 $r_{\rm S}$ for an increasing number of RLAGN samples (Kim et al 2020; Janssen et al 2021; Issaoun et al 2022; Jorstad et al 2023; Fuentes et al 2023; Paraschos et al 2024), covering a variety of jet viewing angles, accretion rates, jet powers and nuclear environments. By comparing the observed jet properties of these sources with those of M 87 at matched gravitational scales, we will be able to better answer the diversity and universality of RLAGN physics across diverse accretion rates. For instance, an ongoing debate is what physics ultimately controls the extension of jet ACZ. As reviewed in Sect. 5, the jet profile of M 87 is constrained very well and many other AGN jets appear to exhibit similar parabolic-to-conical transitions regardless of AGN types, suggesting that such a jet transition is universal. However, while the M 87 jet breaks at the Bondi radius, a majority of other jets seem to break at distances different from their Bondi scales (e.g., Kovalev et al 2020; Boccardi et al 2021; Okino et al 2022). Thus, it is worth continuing to explore how magnetically-driven jets evolve under the stratified atmosphere governed by SMBHs in future observations. This would require more samples with a broad range of AGN parameter space.