© The Authors 2023

1. Introduction

There is increasing evidence that all galaxies with a bulge component have a nuclear supermassive black hole (SMBH; e.g., Saglia et al. 2016). Direct measurements of SMBH masses – available for only ∼200−250 galaxies – are enabled via molecular gas kinematics (e.g., Barth et al. 2016), ionized gas kinematics (e.g., Schnorr Müller et al. 2011), water vapor maser kinematics (e.g., Gao et al. 2017), stellar dynamics (e.g., Rusli et al. 2011), and reverberation mapping (e.g., Bentz et al. 2010). A review of these methods and implications for black hole and galaxy coevolution is available in, for example, Kormendy & Ho (2013).

For this relatively small and biased sample of SMBH measurements, the SMBH mass (M_•) appears to be related to several properties of the host galaxy, such as the bulge luminosity (Kormendy & Gebhardt 2001; Marconi & Hunt 2003; Gültekin et al. 2009) and the bulge stellar velocity dispersion (Ferrarese & Merritt 2000; Gültekin et al. 2009; Saglia et al. 2016). Additionally, the so-called “single-epoch reverberation mapping” method, reverberation mapping of the broad emission line and the luminosity of the nearby continuum in broad line active galactic nuclei (AGNs), can be used to estimate the SMBH mass. These scaling relationships have been used to understand the role of black holes in the formation and growth of their host galaxies (e.g., Di Matteo et al. 2005).

Measurements of a SMBH mass using different techniques can result in inconsistent results (e.g., Walsh et al. 2013; Smith et al. 2021). A famous example is NGC 4486 (M 87), for which the measurement from stellar dynamics (Gebhardt et al. 2011, hereafter G11) is twice the value from ionized gas kinematics (Walsh et al. 2013, hereafter W13), with the two values differing by more than the 3σ of each quoted measurement.

NGC 4486 (M 87) is a giant elliptical galaxy located in the center of the Virgo cluster. It has a prominent relativistic jet that has been well studied on scales of 10 gravitational radii (Event Horizon Telescope Collaboration 2019) out to ∼40 kpc from the nucleus (e.g., Owen et al. 2000). The jet is projected on the sky with a position angle (PA) of 288° (Walker et al. 2018). The orientation and inclination of the nuclear ionized disk was constrained to PA = 45° and inclination (i) ∼ 42° by Walsh et al. (2013). Given this, Jeter et al. (2019) argue that the jet has an inclination of 18° to the line of sight (LOS) and find that the jet and gas disk axes are misaligned by at least 11° (with a typical misalignment of 27°).

There are several reasons to use M 87 as a test case for black hole mass measurement (in)consistencies. First, it is one of the brightest galaxies in the sky. Second, as a massive cD elliptical with a large velocity dispersion, it is expected to host a massive, ≥10⁹ M_⊙, SMBH. Finally, it hosts one of the most studied and well-known “radio-loud” AGNs.

The first attempts to measure M_• in M 87 were made by Sargent et al. (1978) and Young et al. (1978), who suggested a black hole mass of about 10⁹ M_⊙. In companion papers, Ford et al. (1994) presented Hα+[N II] λλ6548,6583 emission-line images of the M 87 nucleus using data from the WFPC2 camera aboard the Hubble Space Telescope (HST), and Harms et al. (1994) presented spectra at six positions in the M 87 nucleus using data from the Faint Object Spectrograph aboard HST. While the nuclear emission line region was clearly not an isolated ionized disk, their isophotal fits to the emission line flux distribution in the nuclear arcsecond, the presence of trailing spiral arms, and the large emission line velocity gradients between apertures led them to posit a nuclear ionized gas disk at PA ∼ 1°–9°, inclination ∼42°, and a black hole mass of 2.4 × 10⁹ M_⊙ assuming a distance of 15 Mpc to M 87. Macchetto et al. (1997) used three long-slit spectra from the Faint Object Camera aboard HST and found that the kinematics are best explained by a 3.2 × 10⁹ M_⊙ black hole (for a distance of 15 Mpc and i = 51°) but noted that the disk inclination is the most uncertain parameter, with likely values between 47° and 65°.

G11 used Schwarzschild modeling of the observed stellar brightness and LOS velocities to measure a black hole mass of (6.6 ± 0.4) × 10⁹ M_⊙ for a distance of 17.9 Mpc. They used data from the Integral Field Spectrograph on the Gemini North Telescope with adaptive optics (AO) correction, combined with extensive kinematics out to large radii. The large-scale stellar kinematics came from two sources: SAURON integral field unit (IFU) data from Emsellem et al. (2004) and VIRUS-P IFU data from Murphy et al. (2011).

Walsh et al. (2013) measured the black hole mass of M 87 via the kinematics of the Hα+[N II] λλ6548,6583 emission lines using data from the Space Telescope Imaging Spectrograph (STIS) aboard HST. Their spectral data were taken in five parallel slits at PA = 51°, separated by 0$\stackrel{″}{.}$1 in the plane of the sky, and with the central slit crossing the galaxy nucleus. They interpreted the ionized gas emission lines as coming from a rotating disk, with PA = 45° and i = 42°. With these parameters, and eliminating data from pixels very close to the nucleus where the Hα and [N II] lines are highly blended, they derived a M_• of 3.5${}_{-0.7}^{+0.9}\times {10}^9{M}_{\odot }$ for a galaxy distance of 17.9 Mpc.

Recently, the black hole mass of M 87 was constrained using data from the 2017 campaign of the (EHTC; Event Horizon Telescope Collaboration 2019). A SMBH enveloped in radio-emitting plasma is expected to have a bright gravitationally lensed ring of emission surrounding the black hole “shadow”. The diameter of the ring is predicted to be ∼9.6 − 10.4 GM_•/c², relatively independent of the SMBH spin (Johannsen & Psaltis 2010; Gralla et al. 2020). The ring diameter measured by the EHTC, together with standard General Relativity, implies a M_• of 6.5 ± 0.2 (statistic) ± 0.7 (systematic) × 10⁹ M_⊙ for a distance to the galaxy of 16.8${}_{-0.66}^{+0.75}$ Mpc. Scaled to this distance, the G11 and W13 SMBH masses are, respectively, 6.14${}_{-0.62}^{+1.07}\times {10}^9{M}_{\odot }$ and 3.45${}_{-0.26}^{+0.85}\times {10}^9{M}_{\odot }$ (Event Horizon Telescope Collaboration 2019). The EHTC result is thus compatible with that of G11 but not W13. More recently, the EHTC used the size of the ring around the SMBH in the Galaxy (Sgr A^*), together with standard General Relativity, and found an excellent agreement between the SMBH mass derived from the black hole shadow and the previous high precision (resolved) stellar dynamics mass (Event Horizon Telescope Collaboration 2022). This result lends additional confidence to the EHTC-derived SMBH mass of M 87.

Most recently, Liepold et al. (2023) studied the stellar kinematics of M 87 with data from the integral field spectroscopy at the Keck II Telescope, concluding that the galaxy is strongly triaxial (there is a misalignment of 40° between the kinematic axis and the photometric major axis from ∼5 kpc outward). They derived a SMBH mass of 5.37${}_{-0.25}^{+0.37}\times {10}^9{M}_{\odot }$ (for a galaxy distance of 16.8 Mpc) using Schwarzschild orbit modeling with a varying mass-to-light ratio (M/L).

To summarize previous SMBH mass measurements (and assuming a distance of 16.8 Mpc to M 87): stellar dynamics and EHTC measurements suggest a more massive SMBH (5.4 − 6.5 × 10⁹ M_⊙), while gas kinematics suggest a significantly lower SMBH mass (3.5 × 10⁹ M_⊙). Two immediate explanations for the discrepancy of the ionized gas kinematic value are a lower inclination for the nuclear ionized disk and a sub-Keplerian rotation of ionized gas (Jeter et al. 2019); both are explored in this work. While we explore a full parameter space of SMBH mass 2 − 8 × 10⁹ M_⊙ and inclination 20 − 50°, some position-velocity (pv) diagrams and residual velocity maps use an illustrative high-mass black hole (HBH; 6.0 × 10⁹ M_⊙), a rough mean value of the three high-mass SMBH measurements, and a low-mass black hole (LBH; 3.5 × 10⁹ M_⊙). Furthermore, we also use illustrative models with disk inclination i = 42° (the value derived in Walsh et al. 2013) and i = 25°. For example, the HBH i25 model is a HBH black hole in a disk with inclination 25°.

A more comprehensive understanding of the complexities of the ionized gas kinematics requires IFU spectroscopy. Indeed, the morphology of the ionized gas in M 87, on scales of arcseconds to tens of arcseconds, is filamentary and complex, with evidence of multiple velocity components. The ionized gas maps of Boselli et al. (2019) show filaments extending from the nucleus up to ∼3 kpc to the northwest and ∼8 kpc to the southeast. They used data from the IFU Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the Very Large Telescope (VLT) in its wide field mode (WFM) to model the ionized gas kinematics in the inner arcminute, and show that the filaments have perturbed kinematics with velocity differences of ∼700−800 km s⁻¹ but a fairly uniform velocity dispersion of ∼100 km s⁻¹ over the filaments.

To better understand the difference between the M 87 SMBH masses measured using ionized gas, stellar dynamics, and EHTC ring size, we obtained new observations of M 87 with MUSE in its AO narrow field mode (NFM) mode. The data have significantly higher spatial resolution than the WFM data used by Boselli et al. (2019) and both a higher spatial resolution and a better signal-to-noise ratio (S/N) than all previous HST/STIS data. Furthermore, unlike the data from STIS, the MUSE NFM cube fully samples the nuclear 8″ × 8″ region and provides a larger wavelength coverage (and thus additional emission and absorption lines) and can be used to model both the stellar dynamics and ionized gas kinematics. Additionally, previous MUSE WFM observations of M 87 over the central 1′ allow constraints to be set on the larger-scale morphology, kinematics, and dynamics.

In this work we use WFM and NFM data to constrain the nuclear ionized gas morphology and kinematics in order to better understand the mismatch between ionized gas and stellar dynamic black hole mass measurements. In a forthcoming work (Osorno et al., in prep.), we will revisit the stellar dynamics black hole mass measurement using the stellar LOS velocity distributions measured in the same datacubes. The dust and ionized gas structures, and their relationship with the prominent jet, will be explored in Richtler et al. (in prep.).

For clarity, we mention the a priori parameters used in our analysis here. As mentioned above, the EHTC adopted a distance to M 87 of 16.8 Mpc (thus a scale of 0.081 kpc per arcsecond), and we use this distance in our work. G11 and W13 used a distance of 17.9 Mpc to derive their SMBH masses; these masses were scaled to the EHTC distance in Event Horizon Telescope Collaboration (2019), and we use these updated values (listed above). The heliocentric recessional velocity adopted by the NED database for M 87 is 1284 km s⁻¹ (from stars). The radial velocity listed in the SIMBAD database is 1256 km s⁻¹. The best fit of W13 to the ionized gas kinematics in HST/STIS data yielded a radial velocity of 1335 km s⁻¹. Our fit to the stellar absorption lines in an annular radius between 2″ and 24″ in the MUSE WFM cube yields a radial velocity of 1313 km s⁻¹. Our equivalent fit to the full field of view (FOV) of the NFM cube yields a radial velocity of 1301 km s⁻¹. The rotational component of the diffuse ionized gas (i.e., not considering the velocity of the spiral arms) in the inner 1″ gas disk in our NFM cube is best centered using a radial velocity of 1260 km s⁻¹. In this work, unless explicitly mentioned otherwise, we use “zero” radial velocities of 1313 km s⁻¹ and 1260 km s⁻¹ for the stellar and ionized gas kinematics, respectively.

The rest of the paper is organized as follows: observations and data processing are presented in Sect. 2; our programs and modeling procedure in Sect. 3; relevant stellar kinematics and properties in Sect. 4; the ionized gas moment maps in Sect. 5; the nuclear ionized disk geometry in Sect. 6; ionized gas filaments in Sect. 7; ionized gas pv diagrams along major and minor axes in Sect. 8; residual velocity maps in Sect. 9; constraints on black hole mass in Sect. 10; the biconical outflow in Sect. 11; and finally the discussion and conclusions of this study in Sect. 12.

2. Observations and data processing

Observations were made with the MUSE IFU mounted on the VLT (Bacon et al. 2010). MUSE covers a wavelength range of 4650−9350 Å with a wavelength sampling of 1.25 Å per pixel in both WFM and NFM modes. The spectral resolution is R = 3026 at 7000 Å. The WFM mode has a FOV of 60″ × 60″, and a spatial sampling of 0$\stackrel{″}{.}$2 per pixel. The NFM mode has a FOV of 8″ × 8″ with a spatial sampling of 0$\stackrel{″}{.}$0252 per pixel. Since NFM data are obtained in AO mode with sodium laser artificial guide stars, no data are obtained in the wavelength range of 5780−6050 Å.

The WFM mode observations were made on June 28, 2014, during a MUSE Science Verification run. Data were taken at airmass 1.3, with a differential image motion monitor (DIMM) seeing of 0$\stackrel{″}{.}$9. The NFM mode observations were made on February 2, 2020, as part of project ID 0103.B-0581 (PI: Nagar). Data were taken at airmasses between 1.25 and 1.5, with DIMM seeing values of 0$\stackrel{″}{.}$4 to 0$\stackrel{″}{.}$5.

The WFM data were processed with the standard MUSE pipeline (v1.6.1) and the calibrated datacubes were hosted on the ESO Archive (Processed Data Science Portal). A single processed datacube, with a total on-sky exposure time of 1800 s, was downloaded from this archive. The NFM data were processed by the standard MUSE pipeline (v2.8.4) and the calibrated datacubes were hosted on the same ESO Science Portal. Three processed datacubes, each with 2100 s of on-sky exposure, were downloaded from this archive. A fourth observation, which was aborted after 83 s of exposure, was not used. The three NFM datacubes were aligned using the position of the bright nucleus and jet knots (one of the cubes had to be shifted by 2 pixels to the E and 1 pixel to the north) and combined into a single 6300 s exposure datacube. The NFM-AO mode point-spread function (PSF) is best defined as a double Moffat function: one for the core and the other for the wings. Given that M 87 is relatively northern for the VLT and has to be observed at a relatively large airmass, the full width at half maximum of the PSF core is expected to be better than 100 mas (see the MUSE Exposure Time Calculator¹ and MUSE User Manual²). With no stars in the FOV, we were unable to measure the true PSF in our data set. However, three globular clusters in M 87 can be distinguished in a continuum image of the NFM cube after subtracting a model galaxy bulge. The full width at half maximum of the profiles of these globular clusters range between 82 and 93 mas.

The final NFM datacube provides full areal coverage of the inner 8″ × 8″, at a spatial resolution better than HST, a depth significantly better than the previous observations of W13, and an almost complete coverage of the optical wavelength range. In Appendix A we compare the two data sets, with Fig. A.1 highlighting the large improvement the new data set brings.

3. Methods and software

In this section we very briefly describe the process and programs used to obtain the stellar and ionized gas kinematics, and the models used in our analysis. A more detailed description of these can be found in Appendix B.

Most codes used, both publicly available and our own, are in Python, and we extensively used the “astropy” library. We used the version 3.1.0. of the GIST Pipeline software (Bittner et al. 2019), which integrates programs for deriving the stellar and ionized gas kinematics of galaxies from optical data cubes. Stellar kinematics, and stellar population properties, were constrained using the version of the penalized pixel fitting technique (pPXF; Cappellari & Emsellem 2004; Cappellari 2017) included in the GIST Pipeline. Different binning techniques were used to bin spaxels in order to increase the S/N of the spectra, and spaxels that include the jet were excluded from the analysis. Details are provided in Appendix B.

Ionized gas kinematics was derived in two ways. First, using GANDALF (Sarzi et al. 2006; Falcón-Barroso et al. 2006; Bittner et al. 2019) within the GIST Pipeline. Here, in each spectrum, a simultaneous fit to both absorption (from a previous pPXF run) and emission lines, is used to obtain an “emission-line-only” spectrum. Since spatial binning was required to increase the reliability of the absorption line fits, the resulting emission line kinematics are obtained in the same bins. While these GIST results were used to check the results below, we do not show them in the figures of this work. Second, to maintain the full intrinsic resolution of the datacube, for each emission line of interest we fit a single Gaussian (for the emission line), plus a first order polynomial (for the continuum), to obtain the intensity, velocity, and velocity dispersion of the emission line in each spaxel. In the case of emission line complexes with wavelengths close enough to be blended in the datacube, for example Hα+[N II] λλ6548,6583, we simultaneously fit each emission line with a single Gaussian. Velocities and dispersions of the different emission lines were tied, and some line ratios were fixed, in these multiple Gaussian fits (see Appendix B).

We also built pv diagrams to reveal the full detail of the ionized gas kinematics, especially in the presence of multiple components. These are extracted from the data cubes along both straight slits and curved pseudo-slits.

In Sects. 8–10, we compare observed velocity maps and pv diagrams directly to analytic kinematic models and to simulated cubes based on analytic kinematic models. In the “Keplerian” model, we assume that the ionized gas is rotating in a thin disk, with a given PA and inclination, under the influence of both the potential of the nuclear black hole and the galaxy. The rotational velocity of the gas (v_gas) is

$$\begin{array}{c}\hfill v_{\mathrm{gas}}(r)=\sqrt{\frac{G[M_{\bullet }+M_{\mathrm{bulge}}(r)]}{r}},\end{array}$$(1)

where r is the nuclear distance, G the gravitational constant, M_• the black hole mass, and M_bulge is the total mass of the stars within a sphere of radius r centered on the nucleus. For its part, M_bulge is derived with the following equation,

$$\begin{array}{c}\hfill M_{\mathrm{bulge}}(r)={\displaystyle {\int }_0^r}4\pi r^{\prime 2}\mathrm{{\rm Y}}(r^{\prime })\nu (r^{\prime })\mathrm{d}{r}^{\prime },\end{array}$$(2)

where Υ is the radially varying M/L of the stars, ν is the volume luminosity density of the stars, and r′ is the integration variable. The V-band luminosity density of M 87, as a function of radius, is taken from Fig. 1 of Gebhardt & Thomas (2009). The values in this plot were extracted with version 3.5.1 of the Plot Digitizer Online Application³. As the luminosity density profile is not an analytic function, the integral in Eq. (2) was performed numerically. Given the large black hole mass, the galaxy potential is significant only at r ≳ 2″. The radially varying M/L ratio in the V-band is taken from Liepold et al. (2023), who fit the R-band M/L variation determined by Sarzi et al. (2018) with a logistic equation, obtaining, for example, values of V-band M/L of 8.65 and 3.46 in the central and the outer part of M 87, respectively.

The LOS velocity of the ionized gas in a given spaxel or along a slit was obtained by projecting v_gas with the following equation,

$$\begin{array}{c}\hfill v_{\mathrm{los}}=v_{\mathrm{gas}}sinicos\theta ,\end{array}$$(3)

where i is the inclination of the disk with respect to the plane of the sky, and θ is the angle (in the plane of the disk) between the major axis of the disk and the deprojected PA of the slit.

We also compared our observables to the radiatively inefficient accretion flow (RIAF) model, which Jeter et al. (2019) advanced as a potential explanation for the ionized gas derived black hole mass being lower than the stellar derived one in M 87. This model considers the presence of sub-Keplerian rotational in the ionized gas, with support compensated with a radial outflow in the disk. Jeter et al. (2019) parameterized the radial and circular velocities as, respectively, v_r = −αv_kep and v_ϕ = Ωv_kep, where v_kep is the Keplerian velocity (Eq. (1)). While α and Ω in RIAFs can have a range of values, we followed the specific RIAF parameters used by Jeter et al. (2019): $\alpha =\sqrt{0.1}$, $\mathrm{\Omega }=\sqrt{0.7}$, and a velocity dispersion profile $\sigma =\sqrt{f}{v}_{\mathrm{kep}}$, with f = 0.1.

We used the program Kinemetry⁴ (Krajnović et al. 2006) to constrain the best-fit PA and inclination of the ionized gas disk. Kinemetry fits the LOS velocity in concentric ellipses, and provides the best-fit PA, inclination, circular velocity, and noncircular velocity components for each fitted ellipse (i.e., at each radius). Velocity field fitting can be done with the PA and inclination as free or fixed parameters.

We created simulated emission line MUSE cubes using the KinMSpy software (Davis et al. 2013). While this was designed to simulate molecular gas cubes, since the third axis of the cube is in velocity space, it is easily adapted to MUSE datacubes. KinMSpy allows our analytical rotation models to be convolved with the same spatial and spectral broadening and sampling of the observed data set. When creating a simulated emission line cube for lines from the nuclear disk, we used the program skySampler in KinMSpy to replicate the surface brightness of the emission lines; skySampler distributes the ∼2 million simulated disk clouds in order to best fit the intensity image of the emission line. The velocity differences between the observed and simulated velocity maps were compared using several metrics and masks in order to obtain the best-fit simulated cube and, thus, the analytic model.

The “toy” nuclear ionized outflow in Sect. 11 is also simulated with KinMSpy. This is modeled using one or two filled cones with vertices in the nucleus. We define our coordinate system with the x-axis positive to the east, the y-axis positive to the north, the observer along the positive z-axis, and the origin at the nucleus. The cones inclination to the observer, their opening angles, lengths, and their radial outflow velocity, were adjusted – by eye – for a good general agreement with the observed velocity maps.

4. Stellar properties from the WFM datacube

A full analysis of the stellar dynamics and Schwarzschild modeling of the LOS stellar velocity distribution to constrain the black hole mass, is deferred to a future work (Osorno et al., in prep.). Here we only present three stellar-related maps most relevant to our interpretation of the ionized gas kinematics.

Stellar moment maps and weighted mean ages and metallicities of the best-fit stellar templates were constructed with pPXF within the GIST Pipeline (see Appendix B for details). The left panel of Fig. 1 shows the velocity map derived from the WFM cube. A systemic (radial) velocity of 1313 km s⁻¹ was used as zero velocity. As seen in the WFM velocity map, and as earlier noted by Emsellem et al. (2014), there are two rotation patterns: one of ±5 km s⁻¹ in the inner ∼18″ with redshift to the northeast, and the other of ±20 km s⁻¹ (in the opposite direction, i.e., redshift to the southwest) and going out to beyond the FOV of the cube. The velocity dispersion map (not shown) is symmetric with values ∼250 km s⁻¹ in a ring of nuclear radius ∼20″, and steadily increasing to more than 400 km s⁻¹ with decreasing nuclear distance. Given that, as expected, stellar dispersion dominates stellar rotation in the nucleus; therefore, we ignored these ±5 km s⁻¹ or ±20 km s⁻¹ of stellar rotational velocities when modeling the gas kinematics in the next sections.

Fig. 1. Maps of stellar properties derived with the GIST Pipeline run on the MUSE WFM datacube. Left: stellar velocity map; zero velocity corresponds to a radial velocity of 1313 km s−1. Middle: map of the weighted mean stellar age of the best-fit templates. Right: map of the weighted mean metallicity of the best-fit templates. These maps were derived using GIST Pipeline over a wavelength range of 5050−6000 Å, with Voronoi binning used to obtain a S/N ≥ 300 in each bin. All spaxels with significant jet emission were masked in the datacube to avoid confusion.

The middle and right panels of Fig. 1 show the maps of weighted stellar age and metallicity from the WFM cube. While the fitted stellar population gets older when approaching the nucleus, the stellar age and metallicity do not vary significantly out to ∼12″, with values between 13 and 13.5 Gyr and 0.32−0.35, respectively (except in the innermost arcsecond, where the large dispersion, and the nuclear jet emission, potentially confuse the fits). Outside ∼12″, the stellar age decreases rapidly with increasing nuclear distance, reaching values below 10.5 Gyr. The metallicity also decreases with radial distance between 12″ and ∼23″ (values between 0.2 and 0.24) and then slightly increases to values between 0.24 and 0.3. It is worth noting the steep change in the metallicity at radius ∼12″, roughly the outer radius of the counter-rotating core. For a V-band M/L ratio of 4 (as used by W13) and a G11 black hole mass, the enclosed galaxy mass equals the black hole mass at radius $r\sim 6\stackrel{″}{.}75$. Our use of a higher V-band M/L ratio in the nucleus, following Liepold et al. (2023), still does not affect our results in the inner arcsecond of the ionized gas disk, and is mainly relevant in the analysis of the kinematics of ionized filaments at ≳2″.

5. Ionized gas: Moment maps

The unprecedented depth, spatial resolution, and full areal coverage of the MUSE WFM and NFM cubes reveal new complexities in the nuclear ionized gas kinematics. Since the [N II] λ6583 emission line is the brightest in both the NFM and WFM cubes, we used this line extensively in our analysis. However, in the innermost ∼1$\stackrel{″}{.}$5, the emission lines are very broad and have multiple velocity components, so we were unable to clearly separate the Hα and [N II] lines in our Gaussian fits. We thus also relied on the [O I] λ6300 emission line in this area, since this is a relatively isolated emission line and is bright enough to derive reliable kinematics in the inner arcsecond.

5.1. Hα+[N II] WFM moment maps

Ionized gas maps of the central arcminute of M 87 were previously presented by Boselli et al. (2019) using the same WFM data cube used here; for these maps we thus only briefly review relevant known results, and discuss features relevant to the black hole mass determination. Figure 2 (top row) show the intensity, velocity, and dispersion maps of [N II] λ6583 obtained from the WFM data cubes. The velocity maps assume a systemic velocity of 1260 km s⁻¹, and a sigma clipping was applied to remove bad pixels. Several features are worth noting in these moment maps. The two dotted black lines, indicated with green arrows, trace the two primary ionized gas filaments that enter the nucleus (in projection). The filament labeled “W1” is the most prominent nuclear filament and extends beyond the MUSE WFM FOV to the west (Gavazzi et al. 2000; Werner et al. 2010; Boselli et al. 2019). The filament marked “W2” also crosses the nucleus but from the east, and connects to filament W1 about 10″ to the northeast of the nucleus. There is an additional highly redshifted connection from the midpoint of filament W2 to filament W1. Additional shorter filaments and clouds, indicated in the velocity map with purple arrows, show velocities relatively different from that of filament W1. The velocity map shows a ∼6″ diameter region to the north of the nucleus (indicated with the red arrow in the velocity and dispersion maps) that is ∼200 km s⁻¹ blueshifted from systemic. We argue below that this is from radial outflow whose axis is aligned with the jet, which adds additional complexity to the nuclear kinematics. In the innermost arcsecond one can already distinguish the posited disk-like rotating component with redshifts (blueshifts) up to ∼400 km s⁻¹ to the northwest (southeast).

Fig. 2. Moment and residual maps for the [N II] λ6583 line from the MUSE WFM datacube. Top: moment maps of the [N II] λ6583 line as derived from Gaussian fits to the datacube. From left to right are the total intensity, velocity, and dispersion maps. Two ionized gas filaments (dotted black lines) are indicated by green arrows and the labels W1 and W2; these dotted lines are the pseudo-slits along which the pv diagrams of Fig. 3 are extracted. Several ionized gas bubbles are indicated by purple arrows, and the position of the outflow is indicated by a red arrow. Bottom: residual velocity maps of the [N II] λ6583 line in the datacube. These maps, discussed in Sect. 7.1, show the residual velocity field after subtracting a model of a rotating thin nuclear disk (in a SMBH plus stellar potential) for a black hole and nuclear disk with the following parameters (M•, disk inclination): 6.0 × 109 M⊙ and 42° (left); 3.5 × 109 M⊙ and 42° (middle); and 6.0 × 109 M⊙ and 25° (right; our preferred model; see Sect. 7.1); all have a disk PA of 45°.

The dispersion map shows high dispersion to the northeast of the nucleus; we argue below that this is due to mixing rotation in the disk with the blueshifted cone of the radial outflow. The dispersion is high in a larger region to the southwest of the nucleus. We argue below that this is due to mixing rotation of the disk with the innermost part of filament W1, which crosses the nucleus in projection but maintains a blueshifted velocity, that is to say, the filament kinematics does not follow rotation in the potential of the black hole. In Fig. 3 we present pv diagrams for the W1 and W2 filaments; they are discussed in Sect. 7.1.

Fig. 3. Position-velocity diagrams of the [N II] λ6583 line for pseudo-slits along the gas filaments W1 and W2 identified in Fig. 2. The x-axis zero position of each filament corresponds to the black cross in Fig. 2. The three curves are the expected disk rotation velocities with the parameters described in Sect. 7.1 and Fig. 2, with colors as listed in the insert. At offsets beyond 10″, the velocities of both filaments are best fit with the i = 25° model. This figure is discussed in Sect. 7.1 but placed here for easy comparison with Fig. 2.

5.2. Hα+[N II] NFM moment maps

Equivalent NFM moment maps are presented in Fig. 4, with the locations of several pseudo-slits identified. The dotted black lines correspond to eight possible ionized gas filaments: most are easily distinguishable in at least one of the three moment maps. The three pseudo-slits shown with dotted magenta lines do not denote filaments; these pseudo-slits cross regions in the map whose velocity features help us better understand the nuclear kinematics. The nuclear patterns seen in the WFM moment maps (Fig. 2) – the blue region to the north, the posited rotating disk, and the multiple filaments – are now resolved in exquisite detail. Given the rich complexity of the components and kinematics, these NFM moment maps are more extensively discussed in Sects. 7.2 and 12. At this point we note the diffuse and constant-velocity blue- (red-) shifted emission to the north (south) of the posited rotating accretion disk, which has high velocities but low dispersion, the multiple filaments to the north (with velocities apparently lower than the diffuse blue emission), and the high dispersion in the posited accretion disk. This high dispersion is caused, in part, from the mixing of multiple velocity components: the disk and outflow (north), and the disk and a filament crossing the nucleus in projection (south). Figures 5–7 show pv diagrams for some of the filaments; they are discussed in Sect. 7.2.

Fig. 4. [N II] λ6583 emission line moment maps from Gaussian fits to the NFM cube. From left to right are maps of the total intensity, velocity, and velocity dispersion. Eight ionized gas filaments (dotted black lines) and three pseudo-slits (dotted magenta lines) are numbered and indicated with arrows of the same color in all panels.

Fig. 5. Same as Fig. 3 but for four of the eight filaments marked with dotted black curves in the NFM moment maps of Fig. 4. The remaining four filaments are presented in Fig. C.1. These pv diagrams are discussed in Sect. 7.2 but are placed here for easy comparison with Fig. 4.

Fig. 6. Same as Fig. 3 but for pseudo-slits C1 and C2, marked with dotted magenta curves in the moment maps of Fig. 4.

Fig. 7. Extra pv diagrams of circular pseudo-slits in the [N II] λ6548 moment maps (Fig. 4) Left: same as Fig. 3 but for the pseudo-slit C3 (nuclear radius 0$\stackrel{″}{.}$75), which is marked with a dotted magenta line in the left and central panels of Fig. 4. Right: same as the left panel but for a pseudo-slit at a nuclear radius of 0$\stackrel{″}{.}$375. The range of velocities is larger in order to show the full Hα+[N II] line complex.

A detailed view of the inner ionized disk as traced by the [N II] λ6583 total intensity is shown in Fig. 8. Ford et al. (1994) used isophotal fits to the inner arcsecond ionized disk (which we show in this work to be a mix of different components, not only a rotating disk) to derive a PA of 1–9° (significantly different from the kinematic PA of ∼40–42° supported by all kinematic studies that allowed this PA to vary) and an inclination of 42°. Here, a complex filamentary structure is still notable, and we thus concentrate only on the spiral arms. The most plausible explanation for these spiral arms are trailing arms in a rotating disk as noted by Ford et al. (1994). They are thus the most reliable indicator of the true extent of the disk, and thus its inclination. As can be seen in Fig. 8, the outer ends of the spiral arms are visually better enclosed in an i = 25° (or even lower), rather than i = 42°, inclined disk.

Fig. 8. Zoomed-in view of the NFM moment 0 image of the [N II] λ6583 line to better illustrate the spiral arm pattern and disk morphology. For reference, the blue and red ellipses show the expected morphology of a thin disk with inclinations of 42° and 25°, the values favored by W13 and this work, respectively.

5.3. [O I] NFM moment maps

In the innermost arcsecond, the large intrinsic linewidths and the presence of multiple velocity components (especially highly blueshifted components) make the automated Gaussian fits to the Hα+[N II] lines relatively unreliable (though these multiple components are still visually distinguishable in pv diagrams). Among the other emission lines in the NFM cube, our best option is the [O I] λ6300 line as it is relatively isolated and has sufficient S/N.

Figure 9 (top row) shows the moment maps for the [O I] λ6300 line in the NFM data cube; the FOV is smaller here as we only show regions with relatively high S/N. The intensity (moment 0) map shows decreasing flux with increasing nuclear distance but does not clearly delineate an axisymmetric inclined disk with a clearly identifiable inclination. For reference we overplot, in blue and red lines, illustrative ellipses expected for a circular disk with inclination 42° (the value determined by W13) and 25° (the value we favor in this work; see Sects. 6 and 10), respectively. The velocity map shows a clear rotational pattern, but with a twist close to the nucleus, and an amplitude that is not symmetric across the minor axis; the blueshifted emission has smaller amplitudes than the redshifted emission in the radius range 0$\stackrel{″}{.}$1–0$\stackrel{″}{.}$3. In the next sections we argue that this is due to a filament that does not participate in rotation and crosses the disk in projection at a constant blueshifted velocity (which is not as blue as the expectation from rotation). The velocity dispersion is also non-axisymmetric: In the highest S/N inner (≤0$\stackrel{″}{.}$2) disk, the blueshifted side of the disk has a larger velocity dispersion, reaching values of ∼600 km s⁻¹, to the southwest. Spectra here are double-peaked with a red component at velocity similar to the redshifted side of the disk, and a blue rotating component (see Fig. 10). The redshifted component could be attributable to the innermost region of the bi-conical outflow (Sect. 11) or a filament crossing the nucleus in projection (Sect. 7.2). A high dispersion region 0$\stackrel{″}{.}$5 to the west of the nucleus, in an area with a lower S/N, is also seen. Finally, a region to the north has the lowest velocity dispersion observed (∼170 km s⁻¹), which likely traces the intrinsic dispersion of the disk (Fig. 10).

Fig. 9. Same as Fig. 4 but for the [O I] λ6300 line in the MUSE NFM cube. Only a subset of the NFM FOV, in which the S/N of the [O I] line is high enough, is shown. For illustration, the blue and red ellipses in the top panels show the expected morphology of a thin disk with an inclination of 42° and 25°, the inclination values favored by W13 and this work, respectively. In the rightmost panels, a cross marks the kinematic center; the apertures marked with black circles and labeled with numbers 1 to 6 are used in Fig. 10.

Fig. 10. [O I] spectral profiles (solid black lines) in selected apertures of the nuclear ionized disk for the apertures indicated in the right panels of Fig. 9, with their single Gaussian fit overlaid in red. The strong line is [O I] λ6300 and the redshifted weaker line is [O I] λ6364. The systemic velocity (dotted black line) and rotation velocities in models with M• = 6.0 × 109 M⊙ and i = 25° (HBH i25; purple) and M• = 3.5 × 109 M⊙ and i = 42° (LBH i42; cyan) are overplotted. Apertures are labeled by their number and their RA and Dec offset from the kinematic center in arcseconds. Apertures were selected to sample portions of the disk with different velocity dispersions and velocity residuals. Aperture 1 is located in the lowest dispersion area of the disk. Aperture 6 shows two components: one centered on the expectation of rotation and the other, with a potentially larger flux, is redshifted 950 km s−1 from the rotating component and 330 km s−1 from the systemic velocity. Apertures to the north (1, 2, and 3) show smaller dispersions, plus a potential component blueshifted by ∼350 km s−1 that is much weaker in flux than the rotating component.

The asymmetries in the dispersion and residual velocity maps of Fig. 9 merit further details, especially given the filaments that cross the rotating disk in projection in the NFM FOV. The right panels of Fig. 9 (i.e., the dispersion and HBH i25 residual maps) show the locations of six apertures whose [O I] spectra are shown in Fig. 10. Apertures 1 and 3 are in the “spiral arms” seen in the residual velocity maps, while aperture 2 is in a region between them. As seen in these three spectra, the central velocity of the best-fit single Gaussian coincides with the peak velocity of the spectrum, confirming that the velocity map from Gaussian fitting reproduces well the central velocity of the strongest emission component. However, all three apertures show a potential weaker (and well separated) component near ∼6320 Å. This weak blueshifted component is at the expected velocity of our posited blue conical outflow to the northeast (Sect. 11). Aperture 1 is in the region of the lowest velocity dispersion of the nuclear ionized disk. It is cleanly, and at high S/N, fit with a single Gaussian with dispersion ∼170 km s⁻¹. This dispersion likely reflects the true dispersion of the disk in the absence of additional nonrotating velocity components. Apertures 4, 5, and 6 are in regions of high velocity dispersion. Aperture 4 seems to have a single component with high dispersion and a redshifted velocity with respect to the predictions of both HBH i25 and LBH i42 models; it is possible that they are a blend of two components with approximately the same intensity. The double component nature is most clearly seen in the spectrum of aperture 6: here a rotating component is distinguishable from a redshifted component: 950 km s⁻¹ from the rotating component and ∼350 km s⁻¹ from systemic; this velocity offset would be expected if the gas were in the red cone of the nuclear outflow (Sect. 11). Aperture 5 (west of the nucleus) also has high velocity dispersion, but the S/N here makes it impossible to tell if it is formed from one or more components. Overall, apart from small velocity offsets from models, regions of the disk show higher velocity dispersion due to additional components, and should be masked when comparing observed velocities to simulated rotating model cubes. Specifically, we mask the high velocity dispersion regions and/or the spiral arm feature regions in Sect. 9.

6. Ionized disk geometry: PA and inclination from Kinemetry

In this section we determine the best-fit PA and inclination of the ionized gas disk using Kinemetry (see Sect. 3). Three different velocity maps were independently input to Kinemetry: the velocity maps from the single Gaussian fitting of [O I] and of Hα+[N II], and also a directly “collapsed” [O I] velocity map. The directly collapsed velocity map is a moment 1 map of the [O I] λ6300 line, created with the spectral-cube package in python, which only included a limited wavelength (thus velocity range); to avoid contamination from the weaker [O I] λ6364 line (63 Å or 3000 km s⁻¹ redward in rest frame), at each spaxel we only included a range of ±40 Å (or ∼±1900 km s⁻¹) centered on the wavelength predicted by the HBH i25 rotational model for that spaxel.

In the first run, for each velocity map, both PA and inclination were free parameters: the best-fit results are shown in the top two panels of Fig. 11. In the case of PA, all three input velocity maps give relatively similar results: the PA starts at value between 58° and 64° at 0$\stackrel{″}{.}$1, decreases out to radius 0$\stackrel{″}{.}$4, and increases beyond this radius. Our kinematic analysis further below concentrates on radii between 0$\stackrel{″}{.}$1 and 0$\stackrel{″}{.}$8: in the inner 0$\stackrel{″}{.}$1 the lines become broad and blended and velocity isophotes are twisted; farther than 0$\stackrel{″}{.}$8 the S/N of disk ionized gas is low and other components (outflow and filaments) dominate. While only the collapsed [O I] velocity field shows a PA more or less constant over a large part of this useful radius range, the mean PA over this radius range is close to 45° for all three velocity fields. The top right panel of Fig. 11 show the equivalent results for the best-fit inclination. None of the velocity fields give a constant inclination with radius, and the inclination evolution appears to be relatively chaotic with values ∼0−45°. The mean inclination values between 0$\stackrel{″}{.}$1 and 0$\stackrel{″}{.}$8 shown are between 17° and 24° for the three velocity maps.

Fig. 11. Results from Kinemetry for the subarcsecond ionized gas velocity maps. The top two panels show results when both disk PA and inclination are allowed to vary, while the bottom two panels show results when the PA is fixed at 45° and only inclination is allowed to vary. Specifically, the top panels show the variation in the best-fit PA (left) and inclination (right) with nuclear radius. The bottom-left panel shows the variation of the inclination when the PA is fixed to 45°. The velocity maps used were the single Gaussian fit velocity maps of [O I] (red) and Hα+[N II] (green) and the “collapsed” moment 1 map of the [O I] line (black). In each of these three panels, the horizontal lines in the corresponding color show the mean value over radii between between 0$\stackrel{″}{.}$1 and 0$\stackrel{″}{.}$8. The bottom-right panel shows the rotation curve (i.e., the “k1” parameter from Kinemetry; connected blue dots) derived by Kinemetry for the Gaussian-fit [O I] velocity map, and compares this to expectations of Keplerian rotation for the HBH i25 (solid blue line) and LBH i42 (dotted blue line) models, all following the left-hand y-axis. The red line, which follows the right-hand y-axis, shows the k5/k1 ratio (i.e., noncircular velocity to circular velocity ratio) derived by Kinemetry: low values here imply that disk velocities are dominated by rotation. At very small radii, too few pixels are included in each ring and Kinemetry results are unreliable.

Given the results of the first run of Kinemetry (and the value previously derived by W13), we decided to set a fixed PA of 45° and leave only the inclination as a free parameter in the second run of Kinemetry. The bottom-left panel of Fig. 11 shows the result of these runs. The variation in inclination with radius is now slightly less chaotic than in the first run but still varies between 20° and 45°, when it is not 0°. The mean values over the 0$\stackrel{″}{.}$1–0$\stackrel{″}{.}$8 range are remain between 17° and 24°, except for the collapsed [O I] map, which prefers a higher mean value of 33°.

While these Kinemetry results do not concretely point to a disk inclination close to 25°, they provide an additional indication that inclination values closer to 20°–25° are more plausible than 42°. It is also an indication of the relatively complex gas motions in the disk could preclude an accurate measurement of the black hole mass.

7. Ionized gas filaments

In this subsection we present and discuss the kinematics of the ionized gas filaments seen in the outer (WFM) and nuclear (NFM) regions. We present evidence that each has its own flow velocity, and only some are affected by the SMBH potential, while others cross the nucleus only in projection. That is, the nuclear “gas disk” is not an isolated relaxed and well-ordered rotating disk, but a mix of gas rotating in a “disk”, emission from gas in filaments, some of which enter with their own proper bulk velocity, and some of which only cross the nucleus in projection. Further below, we also discuss the extra complexity of disentangling the nuclear biconical outflow.

7.1. Wide field mode filaments

In Sect. 5.1 we highlighted two large-scale filaments in the WFM moment maps (Fig. 2). The pv diagrams of the [N II] emission line along these filaments are shown in Fig. 3, together with the predictions of disk rotation in the potential of the SMBH and the galaxy under different assumptions of black hole mass and disk inclination. Since the apertures are curved, we refer to these as pseudo-slits. The offsets (x-axis) in the pv diagrams start at the position of the nucleus, and continue along the dashed lines shown in Fig. 2. That is, the x-axis value indicates the distance traveled along the line from the nucleus (indicated with a black cross in Fig. 3) and not the radial distance from the nucleus. The colored solid lines in the pv diagrams denote the velocities expected from gas in a rotating thin disk for the following models: LBH i42 (magenta), HBH i42 (white), and HBH i25 (green).

All of these models have a disk at PA = 45°. The galaxy potential is obtained as described in Sect. 3, and becomes significant only at radii beyond ∼2″. The LBH i42 and HBH i42 (magenta and white) models use the inclination and PA found by W13 for the ionized gas disk. The green line is illustrative for a HBH and a disk with the lower inclination suggested by our data (see Sect. 8).

Filament W1 is the longest and most prominent. None of our velocity models follow the observed velocities in the first 14″. The filament leaves the nucleus at the PA roughly close to the minor axis of the W13-posited gas disk and one would expect relatively small Keplerian velocities. Nevertheless, velocities are between ±100 km s⁻¹ (with an extreme value of ∼400 km s⁻¹) over the first ∼30″ of the filament. Between 20″ and 30″ there are two – potentially bubble-shaped – clouds at 90 km s⁻¹ and −350 km s⁻¹. The latter has a velocity highly different from the rest of the filament and our models; this bubble is either not part of the filament (but in the same LOS) or is caused by a kink or instability in the filament. From 36″ onward, the observed velocities follow the overall pattern expected from the HBH i25 rotating disk model in both values and shape. To make the HBH i42 and LBH i42 models agree with the observed velocities over this range, one would have to change the ionized gas zero velocity by ∼80 km s⁻¹ to the blue, which would imply an even greater difference between the recessional velocities of ionized gas and stars (∼1300 km s⁻¹).

Filament W2 is shorter and shows relatively smooth velocity changes along its length, with a relatively low dispersion (compared to filament W1) at offsets larger than 3″. Similar to filament W1, the observed velocities close to the nucleus are systematically bluer than all our rotating disk models, especially in the first 5″. Once more, at distances larger than 10″ the HBH i25 model most closely fits the data in values and shape. It is also interesting to note that filament 2 appears to (at least in projection) connect to filament W1 through an additional highly redshifted bridge to the east of the nucleus (Fig. 2), though we find no explanation for this velocity pattern.

Except in the innermost ∼5″–12″, the observed pv diagrams of both filaments W1 and W2 are well explained – in velocity and velocity gradients – by circular rotation of gas in a disk with inclination close to 25° (green line in the pv diagrams). This can also be seen in the residual velocity maps (bottom row of Fig. 2). At these large radii, the galaxy potential, rather than SMBH mass, drives the circular velocity. This is clearly illustrated by the fact that the HBH i42 (white lines) and LBH i42 (magenta lines) models result in very similar velocity predictions. The disk inclination is thus the dominant parameter, and the degeneracy between black hole mass and inclination when fitting the kinematics of the inner arcsecond of the ionized gas disk no longer plagues us.

7.2. Narrow field mode filaments

We now focus on the eight filaments and three pseudo-slits identified on the higher resolution NFM moment maps (Fig. 4). The pv diagrams of the [N II] λ6583 line along these filaments and pseudo-slits are shown in Figs. 5, 6, and C.1. We recall that the x-axis of these figures are offsets along a curved slit, and not the radial distance from the nucleus. Most of the filaments are to the north of the nucleus, and cross the highly blueshifted region. The pseudo-slits (magenta lines in Fig. 4) do not follow gas filaments, but their pv diagrams are useful to understand the complex nuclear kinematics, so they are also presented here.

Each of the eight filaments in Fig. 5 can be differentiated from the others in at least one of the three moment maps in Fig. 4. We expect to see significant effects of the black hole potential on the ionized gas as we get closer to the nucleus, especially in the inner few arcseconds. However, the filament pv diagrams (Figs. 5 and C.1) are almost all significantly different from all of the model predictions of thin disk rotation in the black hole and galaxy potential. Only filament 8, which is to the west of the highly blueshifted region (see the middle panel of Fig. 4) approximately follows our rotating disk model predictions.

Filament 1 can be distinguished in the intensity map; its pv diagram (Fig. C.1) shows differences of up to ∼400 km s⁻¹ blueward from the rotation model predictions. Filaments 2 and 3 have a semicircular form in the moment maps, with their extremes closest to the nucleus. Both can be distinguished in the intensity map, with filament 2 also partially distinguishable in the velocity and velocity dispersion maps. For both filaments, the disagreement between observed velocities and the rotational models increases with increasing nuclear distance. The disagreement is nevertheless smaller close to the nucleus than the case of filament 1. Filaments 4 and 5, clearly distinguishable in the intensity maps, are in the highly blueshifted zone to the northeast of the nucleus. They are the most blueshifted filaments in the inner nucleus, with blueshifts of ∼250−350 km s⁻¹ from systemic, and more than 500 km s⁻¹ offset from the rotation model predictions. Nevertheless, their dispersion velocities are ≤150 km s⁻¹ and both appear to have only a single velocity component. Filaments 6 and 7 are distinguishable in the intensity and velocity maps. With velocities between −200 km s⁻¹ and −150 km s⁻¹, they are less blueshifted than filaments 4 and 5. Filament 6 shows a higher dispersion (∼200 km s⁻¹) than filaments 4, 5 and 7. Filaments 4 to 7 in the NFM can also be distinguished in the WFM maps: at the lower resolution of WFM, they are distinguishable as three filaments (or one filament and a loop).

Filament 8, the westernmost filament, is clearly distinguished in the intensity and velocity map. It is the only filament to approximately follow our rotating disk models, though it is slightly more blueshifted than the models at offsets of $\sim 0\stackrel{″}{.}1-0\stackrel{″}{.}2$ and $\sim 1\stackrel{″}{.}8-2\stackrel{″}{.}2$. Its dispersion is also small (≤150 km s⁻¹). Given its velocity and morphology, it is likely the inner region of filament W1 of the WFM (Fig. 2).

Overall, filaments 1 to 7 follow individual flow velocities, which cannot be explained by a rotating thin disk. Only filament 8 approximately follows the expectations of a rotating thin disk, but since this filament approaches the nucleus in a PA close to the minor axis of the disk, it cannot be used to constrain the black hole mass.

We now discuss the pv diagrams (Fig. 6) of the pseudo-slits C1 to C3 in the NFM cube (dotted magenta curves in Fig. 4). The pseudo-slit C1 (left panel of Fig. 6) cuts through most of the northern filaments from east to west while maintaining a roughly constant radial distance from the nucleus. Nevertheless its pv diagram shows a constant blueshift of ∼200 km s⁻¹ from systemic except over offsets of ∼1″ − 3″ along the pseudo-slit (effectively the northeast region) where a constant blueshift of ∼300 km s⁻¹ is seen. While this is visible in the velocity map, the pv diagram emphasizes the relatively uniform velocity of the northern filaments, the uniformly small dispersion velocity, and the lack of additional velocity components. Effectively the filaments are bright concentrations embedded within, and sharing the velocity of the more diffuse blueshifted emission region. The observed velocities are even further blueshifted from the rotation model predictions, and in fact appear as a rough mirror image of these.

The pseudo-slit C2 starts in the redshifted side of the inner disk, travels to the southeast along a small filament (see the left panel of Fig. 4), and then crosses the redshifted ionized gas to the south in the NFM moment map. In its pv diagram (right panel of Fig. 6), we note the high dispersion in the first arcsecond, and the fact that the kinematics roughly follow the LBH i42 and HBH i25 rotating disk models in the inner two arcseconds. We see a small velocity discontinuity at ∼2$\stackrel{″}{.}$5–3″ offset along the pseudo-slit, which delineates the separation point between the rotating (inner) and nonrotating (outer) components of the gas. However, the continuity of the velocity transition from one component to the other indicates interaction between the two components.

The pseudo-slit C3 is a circle with radius 0$\stackrel{″}{.}$75 centered on the nucleus. The x-axis starts (offset = 0) west of the nucleus and moves clockwise. Over all of the pseudo-slit except the northwest quadrant the observed velocities are well fit with either the LBH i42 model or a HBH i25 model. In the northwest quadrant, the observed velocities are systematically ∼200 km s⁻¹ blueshifted from the model predictions. While filaments 6 and 7 of the NFM can be seen as separate components at even larger blueshifts, the ∼200 km s⁻¹ blueshift in the northwest does not appear to be an additional component but rather a warp or systematic bulk velocity in the disk. At this radius, the southern (especially the southwestern) section of the disk appears the least confused with additional components or velocity offsets. The pv diagram shown in the right panel of Fig. 7 is for a circular slit like C3, but at half its radius (now $r=0\stackrel{″}{.}375$). Here, the disk velocities are more symmetric in all PAs, with the LBH i42 and the HBH i25 both providing a reasonable fit. We note that the northern blueshifted filaments are still clearly distinguishable. At this smaller radius, the gas to the southwest shows a second velocity component: a nonrotating component ∼300 km s⁻¹ blueshifted from systemic velocity. This component dominates the rotating component at smaller radii, and is the primary reason why the blueshifted side of the rotating disk shows smaller absolute velocities than the redshifted side of the disk over radii 0$\stackrel{″}{.}$1–0$\stackrel{″}{.}$3. The velocity of this nonrotating component matches that of filaments 6 and 7 as they approach the nucleus in projection.

8. Ionized gas rotation: Position-velocity diagrams and comparison to analytic models

Figure 12 shows the pv diagrams of the Hα+[N II] λλ6548,6583 line along nuclear slits at PA = 45° (the W13-based disk major axis) and 135° (disk minor axis) in the WFM cube. Even at this relatively low spatial resolution, and at large nuclear offsets, the complex kinematics and presence of multiple components is already obvious. The inner 3″–4″ to the northeast shows highly blueshifted (instead of redshifted) emission, which, with increasing radius smoothly transitions to systemic velocity. Furthermore, at 8″–16″, the bright emission line gas does not follow rotation in a disk with 42° inclination. Along the minor axis to the northwest, the nuclear gas is also blueshifted, and other kinematic “wiggles” are seen further away from the nucleus. Figure 13 shows the equivalent pv diagrams in the NFM cube (left and middle panels) together with the pv diagram of the [O I] λ6300 line along a nuclear slit at PA = 45°. In the left panel, the gas ∼1$\stackrel{″}{.}$5–4″ to the northeast of the nucleus (positive offsets in the figure) primarily traces the kinematics of the filaments and outflow, and there is no clear sign of rotating gas at these offsets. When the nuclear offset decreases below 1$\stackrel{″}{.}$5, instead of continuing “into” (in projection) the nucleus at the same bulk velocity offset, the gas velocity gets redder very rapidly but smoothly (crossing systemic velocity at 0$\stackrel{″}{.}$8 in projection) and then further increases until it is indistinguishable from gas in the rotating disk. Adding a simple bulk (systemic) blueshift to our rotation models cannot reproduce this velocity pattern (the observed velocities drop too rapidly with increasing nuclear distance). A filament crossing the nucleus only in projection would also not explain this pattern, as the projected nuclear separation is always larger than the true three-dimensional nuclear distance. For illustration the dotted green line shows the predictions of gas rotating around a G11 black hole, but at an inclination of 60° and with a bulk blueshift of 650 km s⁻¹. A possible explanation for this feature is that nuclear gas, initially feeding the black hole in a plane with a larger inclination than the main disk, is entrained into the blue cone of the conical outflow (see Sect. 11).

Fig. 12. Position-velocity diagrams of the Hα+[N II] λλ6548,6583 lines in the WFM cube along the W13 major (45°, left) and minor axes (135°, right). The models (see the inset) and the zero velocity are plotted with respect to the [N II] λ6583 line.

Comparing the observed pv data to the models, the HBH i42 overpredicts the observed velocities, while both the LBH i42 and the HBH i25 models are relatively indistinguishable when compared to the data. Along the minor axis (central panel of Fig. 13), the blueshifted gas to the northwest is clearly visible at an almost constant ∼200 km s⁻¹ blueshifted velocity. The [O I] pv diagram (right panel of Fig. 13) covers a smaller range of offsets. Here the HBH i42 model is clearly in the upper envelope of the observed velocities.

Fig. 13. Position-velocity diagrams of the Hα+[N II] λλ6548,6583 (left and central) and [O I] (right) emission lines in the NFM cube, along the W13 major (45°, left and right) and minor axes (135°, central). Rotation models overplotted follow the inset. The dotted green line is for an illustrative rotating model with M• = 6.0 × 109 M⊙ (HBH), an inclination of 60°, and a constant blueshift of 650 km s−1 from the systemic velocity.

All WFM and NFM pv diagrams clearly show filaments nearing the nucleus (in projection) but not following rotational kinematics of the black hole. In several cases, multiple velocity components are seen, so that single Gaussian fits are not usable here. While the highly blueshifted components can be disentangled in pv diagrams at most PAs, for reasons of brevity we cannot show all these here. To avoid some of the confusion created by the highly blueshifted filaments, we instead used the [O I] λ6300 emission line, which, while weaker, is relatively isolated from other emission and absorption lines and is detected out to a ∼1″ radius (right panel of Fig. 13).

9. Residual velocity maps

The results of our Kinemetry (Sect. 6) and residual (observed minus simulated) velocity map analysis (Sect. 10) support low inclinations for the ionized gas disk. To better visualize the differences when using low- and high-inclination disks, and to discern residual velocity features after rotation velocity subtraction, this section presents the residual (observed minus model) velocity maps of competing rotational models. The model velocity maps used in this section are created using KinMSpy.

It is important to note that the RIAF model includes a radial outflow in the disk plane, which results in an observed kinematic PA that is offset from the rotational kinematic PA. For this specific RIAF used here, a rotation with a major axis at PA = 27° added to the outflow component results in an overall observed kinematic PA of 45°. Alternatively, if the near and far side of the disk were exchanged, then a rotation PA of 63° would give the same result.

The bottom panels of Fig. 2 show the residual velocity maps of the [N II] λ6583 line after subtracting the LBH i42, HBH i42, and HBH i25 models. At the relatively large nuclear distances seen in the WFM, the galaxy potential dominates the black hole potential, thus avoiding the degeneracy between inclination and black hole mass in the inner arcsecond. Both the HBH i42 and LBH i42 models (i.e., both models with i = 42°) over-subtract the velocities seen over most of the filaments W1 and W2. The i = 25° model however produces a velocity residual closer to zero. Here it is important to note that we used a recession velocity of 1260 km s⁻¹, which is the velocity that best fits the inner ionized gas disk, and which is 75 km s⁻¹ bluer than the value used by W13. With our recessional velocity of 1260 km s⁻¹, the observed velocity fields of filament W2 (and a large part of filament W1) of the WFM are primarily at velocities between systemic and 180 km s⁻¹ (redshifted). As discussed in Sect. 7.1, reducing the velocity residuals of the i = 42° models requires using a smaller recessional velocity, which would be even further away from the stellar recessional velocity, and thus more unlikely to be correct. In any case, except for some prominent kinks and bubbles in filament W1, and extreme redshifts in the apparent bridge between the W1 and W2 filaments, both W1 and W2 appear to be primarily composed of gas in a rotating (and partially filled) ionized gas disk. The large blue region 1″–10″ north of the nucleus, and the smaller red region ∼2″ south of the nucleus, seen in the residual maps can be explained by the bi-conical outflow (Sect. 11).

Residual velocity maps for the [O I] line are shown in the bottom row of Fig. 9: from left to right are the residuals after subtracting the HBH i42, LBH i42, and HBH i25 models. The HBH i42 significantly overpredicts the LOS observed velocity field. The LBH i42 and HBH i25 models give roughly similar residuals, though the latter provides the velocity residuals closest to zero.

The twisted ionized gas velocity isophotes in the innermost 0$\stackrel{″}{.}$1 (most clearly discerned in Fig. 9, but also clear on close inspection of Fig. 4) is an important result of this work. At a radius of 0$\stackrel{″}{.}$3, the PA of the disk appears to be ∼45° as found by W13 and Kinemetry (see Sect. 6). By ∼0$\stackrel{″}{.}$1 the PA appears closer to 90°. In the innermost pixels the PA appears to twist even more. This innermost twist is best visualized in the residual velocity map (bottom-right panel of Fig. 9). Here the residual map implies a rotation major axis of ∼18°, which is perpendicular to the PA of the inner jet of M 87, as would be expected if the jet was launched perpendicular to the accretion disk. Alternatively, this relatively north-south inner velocity gradient could trace the bases of the bi-conical outflow discussed in the Sect. 11. Macchetto et al. (1997) previously noted the velocity field twist inside 0$\stackrel{″}{.}$2 when modeling their HST spectra with thin disk models. The twist in the apparent PA appears to continue outside the radius of 0$\stackrel{″}{.}$3, as supported by the results in Sect. 6. However, at even larger radii, as shown earlier in this section, the velocity field of the filaments in the WFM maps at r ∼ 10 s of arcseconds can still be fit well with a rotating disk at PA ∼ 45° and inclination 25°.

Emsellem et al. (2014) have shown that the stellar velocity field in the central ∼15″ of M 87 is twisted, with the stellar kinematic PA changing rapidly from 17° to −124°. In the innermost ∼5″, the stars rotate with a major axis close to PA = 20°, with the northeastern side redshifted. This pattern can also be seen in our WFM stellar map (Fig. 2), where, visually, the inner ∼10″ stellar kinematics support a PA of 30° −45° with the northeast receding, which is similar to the pattern seen in the ionized gas here. The main difference is that the stellar rotation velocity amplitudes are relatively low (≲20 km s⁻¹) and that the twisted isophotes are seen on much larger scales.

The other important feature to note in the HBH i25 residual velocity map of Fig. 9 is a three or four armed spiral pattern at ∼100 km s⁻¹ redshift, embedded in a blue (outflow or filaments) halo. This velocity-traced spiral arm pattern closely matches the spiral arm pattern noted by Ford et al. (1994) in their HST ionized gas maps, and seen clearly in our moment 0 maps (e.g., Fig. 8). The spiral pattern is reminiscent of the ionized gas spirals around Sgr A*, which extend out to ∼3 pc, and appear to trace Keplerian circular or elliptical orbits around Sgr A* (Zhao et al. 2009).

10. Observed versus KinMSpy models of gas rotation: Constraining black hole mass and inclination

To obtain a best-fit model (in black hole mass and disk inclination) to the nuclear Keplerian rotating ionized gas disk seen in the NFM cubes, we used KinMSpy to construct simulated datacubes of the [O I] line in the inner 1″, for a range of black hole masses (2 × 10⁹ M_⊙ to 8 × 10⁹ M_⊙ in steps of 0.5 × 10⁹ M_⊙) and disk inclinations (20° to 50° in steps of 2°). In all models, the disk PA was fixed at 45°, and the same galaxy potential (see Sect. 3) was used.

As described in Sect. 3, KinMSpy was given the intrinsic spatial and spectral resolution of MUSE, the pixel sampling of the MUSE cubes, and clouds were placed to best reproduce the respective moment 0 image of the emission line, in order to make the most realistic simulated cubes. Velocity fields of these simulated cubes were then created in python using the spectral-cube package. The residual (observed minus simulated) velocity map was then used with several masking schemes and several metrics, to determine the simulated velocity map that best fits the observed velocity map.

As emphasized in several previous sections, both the nuclear ionized disk and the large-scale ionized gas emission show complex kinematics, often due to the superposition of several distinct components. The masking used on the residual velocity field before calculating best-fit metrics is thus important. We used several masks: (a) Gauss-fit mask: all spaxels that show large residuals between the observed spectrum and the best-fit Gaussian were masked. This effectively masked both low S/N spaxels and spaxels with complex profiles that could not be fit with a single Gaussian; (b) sigma-clip mask: all pixels in which the [O I] dispersion map (Fig. 9) has values larger than 350 km s⁻¹ were masked. This masks also effectively avoids the inclusion of spaxels with multiple velocity components; (c) spiral-arm mask: all pixels that fall within the mini spiral arms seen in the [O I] velocity residual map for the HBH i25 model (lower right panel of Fig. 9) were masked. Specifically we masked pixels with absolute residual velocities larger than 100 km s⁻¹; (d) annular mask: all pixels outside the annular region between 0$\stackrel{″}{.}$2 and 0$\stackrel{″}{.}$4 from the nucleus, and pixels inside a rectangular region between −0$\stackrel{″}{.}$6 and 0$\stackrel{″}{.}$1 in east-west offset, and −0$\stackrel{″}{.}$45 and 0$\stackrel{″}{.}$15 in north–south offset, were masked. This mask retains the highest S/N regions of the disk while deleting the innermost region (with isophotal twists) and the high velocity dispersion nuclear region. The sigma-clip mask is shown in Fig. 14; the remaining masks can be found in Appendix D.1.

Fig. 14. Examples of masks applied to the subtracted velocity maps. The masks were made by removing pixels with velocity dispersions larger than a given limit (sigma-clip mask). From left to right are the masks for [O I] λ6300 and Hα+[N II] λλ6548,6583. See Fig. D.1 for other types of masks.

At each value of SMBH mass and inclination, we calculated several metrics from the masked velocity residual (observed minus simulated) map; multiple metrics were used in order to be robust against changes in the histogram profile and average velocity in a given residual map. These included the standard deviation, the mean of the absolute difference, the difference between the 84 and 16 percentiles, and the r squared value (the last constructed using the observed and model velocities rather than the distribution of the residual velocities). These metrics were then evaluated across the parameter space of SMBH mass and inclination. Most combinations of masks and metrics deliver similar results with respect to favored pairs of SMBH mass and disk inclination, so we show only selected examples here. The exercise above was also made for a comparison of the observed [O I] velocity map with simulated cubes of the RIAF model of Jeter et al. (2019, see our Sect. 3).

Figure 15 presents the resulting metric obtained from the [O I] emission line residual velocity maps for a Keplerian rotation model (left panels) and for the Jeter RIAF model (right panels). Specifically, here we show the results of the standard deviation of the residual velocity maps for three different masks: Gauss-fit, sigma-clip, and spiral-arm. We concentrate first on the left panels in Fig. 15 (i.e., Keplerian rotation models). In general the minimum values of standard deviation follow a “banana” pattern, with almost equal support for high black hole masses at low inclinations as for low black hole masses at higher inclinations, with black hole mass decreasing with increasing inclination. The uppermost panel is the result from the Gauss-fit mask. It favors a high black hole mass and small inclination, with a minimum beyond the range of the plot; the masking here does not necessarily exclude all high dispersion regions and spiral arms seen in the moment (and residual) maps. The sigma-clip and spiral-arm masks (central and lower panels) yield similar results. Both plots show a banana-shaped dark blue region of minimum values, where the metric does not vary significantly. For our specific LBH and HBH values, the inclinations at minimum standard deviation are, respectively, ∼37° and ∼28°. There is very slight support for a HBH at a lower inclination: the standard deviations (for the sigma-clip mask) for the LBH with i = 37° and HBH with i = 28° are 62.2 and 64.8, respectively, a difference of 4.2%; for the spiral-arm mask, the difference is 3.92%, in the same sense. In general, over most masks and metrics, a HBH at a low inclination is favored over a LBH at a high inclination only at the ∼4% level.

Fig. 15. Standard deviation (in km s−1) of values in the residual (observed minus KinMSpy model) velocity map, as a function of the black hole mass and disk inclination of the model, following the color bar to the right of each panel. The left column is for the Keplerian disk model, and the right column for the Jeter RIAF model. From top to bottom are the standard deviations measured using Gauss-fit, sigma-clip, and spiral-arm masks. The panels of each row follow the same color bar.

The equivalent results for the RIAF model residuals are shown in the right panels of Fig. 15. As with the Keplerian model, we see a dark blue banana-shaped minimum region, but this is offset with respect to the Keplerian results as expected: sub-Keplerian rotation results in higher inclinations favored for a given black hole mass. Nevertheless, the level of sub-Keplerian rotation in this specific RIAF model is not sufficient to support the scenario of a HBH at inclination 42°: values of Ω significantly lower than $\sqrt{0.7}$ would be required. Once more the Gauss-fit mask (uppermost panel) favors a high black hole mass and low inclination beyond the plot. For the sigma-clip and spiral-arm masks the inclination ranges for a minimum standard deviation for the LBH and HBH are, respectively, ∼40°–43° and ∼30°–31°. The respective standard deviations are 76.5 and 68.7 for the sigma-clip mask, and 56.4 and 54.4 for the spiral-arm mask. The standard deviation for the HBH case is thus lower by 11% for the sigma-clip mask and 4% for the spiral-arm mask. Thus overall, within the Jeter RIAF model, a HBH in a low-inclination disk provides a slightly better fit than an LBH in a higher-inclination disk, and a HBH in a disk with inclination i = 42° is not supported.

Inter-comparing Keplerian disk scenarios with the Jeter RIAF scenario, we note the following: (a) for a HBH, the Keplerian disk scenario gives lower velocity residuals than the RIAF model across all masks and metrics; (b) interestingly, a LBH in a i = 42° disk is better supported in a RIAF scenario as compared to the Keplerian disk scenario, and larger values of Ω would make the RIAF model even more favorable.

We repeated the exercise for the velocity map of the large-scale filaments seen in the WFM cube, and the masking (a) and (b) explained for [O I]. Here we used the velocity map of the Hα+[N II] line created via Gaussian fitting (with a single Gaussian fit to each of the Hα and [N II] lines). The only differences from the process above is that (a) we varied the inclination between 0° and 50°; (b) the sigma-clip mask was created by masking pixels with velocity dispersion above 100 km s⁻¹; and (c) we only used the Keplerian rotation models. The resulting standard deviation in the residual velocity map, for intensity masking (left) and velocity dispersion masking (right) are shown in Fig. 16. In both cases, the minima are seen at relatively small inclinations, almost independent of the black hole mass. This is since most of the non-masked area is outside of the sphere of influence of the SMBH. The minimum standard deviation is seen at an inclination of about 18° for the Gauss-fit mask and 24° for the sigma-clip mask; using the other metrics and masks give similar results. This supports the scenario that the filaments are rotating in a disk with inclination ∼20°, close to the value we support for the inner subarcsecond disk.

Fig. 16. Same as Fig. 15 but for the WFM Hα+[N II] λλ6548,6583 velocity map (i.e., the extended filaments) and the Keplerian rotation model. Both panels show the standard deviation of the residual velocity map, with the left panel showing the result with the Gauss-fit mask and the right panel the sigma-clip mask.

11. Ionized gas outflow

The moment maps and pv diagrams in the diffuse blue region and filaments to the north of the nucleus clearly demonstrate the presence of a nonrotating component. Additionally, some excessively redshifted gas is seen ∼3″ southwest of the nucleus. With projected velocities up to 400 km s⁻¹, the most viable explanation for these velocities is an AGN-driven outflow.

The left panel of Fig. 17 shows the residual velocity map (observed velocity minus our HBH i25 rotating disk model) of the [N II] λ6583 line, and emphasizes better not only the blue and redshifted regions to the north and south, respectively, but also the velocity gradients seen in these.

Fig. 17. Modeling of the ionized gas outflow in M 87. Left panel: residual (after subtracting the HBH i25 model for the nuclear disk) [N II] velocity map. We emphasize the gradients in the blue- and redshifted velocities to the northeast and south. Middle and right panels: models of a filled conical outflow made with KinMSpy, with each cone shown separately. The cone half-opening angles are 45° (approaching; middle panel) and 30° (receding; right panel). Each cone has an extent of 3″ along its axis, and the radial outflow velocity within each cone is 400 km s−1. The inclination of both cone axes to the LOS is 18°. The green line shows the projection of the cone axes on the plane of the sky; these project to PA 288°, the PA of the jet axis. The dotted orange lines delineate the projected shape of a disk perpendicular to the cone axes.

One could attempt to reproduce these velocity features with a (fully filled) biconical outflow. Illustrative examples of these are shown in Fig. E.1, where the cone axes were chosen to illustratively “match” the blueshifted and redshifted regions, or to project onto the PA of the jet. Neither of these illustrative models can satisfy the observed velocities while maintaining an origin in an AGN-produced outflow along the axis of the jet.

Instead we invoke a biconical outflow in which the cones are only partially filled with ionized gas; a scenario supported by the filamentary and patchy nature of the ionized gas in the inner 30″. To illustrate this scenario, Fig. 17 shows the velocity structure of each cone individually. Here both cones have the same opening angle (45°), and their orientation is chosen so that both cone axes project to PA −72° (the PA of the jet). The velocity gradients in each bicone match those seen in the residual velocity map (left panel of Fig. 17). Thus, supported by the non-axisymmetric filamentary and patchy ionized gas seen on scales larger than 1″ in both the NFM and WFM cubes, we posit that the diffuse and filamentary ionized gas to the north (outside the accretion disk and within the NFM FOV) is primarily between the nucleus and the observer, and the ionized gas to the southwest of the nucleus (outside the accretion disk and within the NFM FOV) is behind the nucleus from the observers point of view. This partial filling of the approaching and receding cones then well recreates the velocity residual map of Fig. 17. This posited partially filled biconical outflow could also explain two other ionized gas blobs to the south and southeast beyond the NFM FOV, but observed in the WFM FOV residual velocity map (bottom-right panel of Fig. 2).

12. Discussion and conclusions

Our full areal coverage of the nuclear (new deep MUSE-AO NFM cube) and larger-scale (archival MUSE WFM cube) ionized gas allows improved constraints to be set for the morphology and kinematics of the ionized gas in the inner 0$\stackrel{″}{.}$1 to 30″ in M 87. On the larger scales, filaments W1 and W2 (which are tens of arcseconds in extent) show velocities that can, in most of their extent, be fit with a partially filled gas disk at inclination i ∼ 25° rotating in the potential of the galaxy. Nevertheless, certain sections of these filaments show large kinks and highly dispersive velocities. The western arm of W1 appears to cross the nucleus in projection and at a relatively low rotation-corrected velocity. While it cannot be clearly traced between 1″ and 2″ west of the nucleus, its incoming direction aligns with a high velocity residual region (after rotation correction) in the southwestern quadrant of the arcsecond-scale ionized disk. We speculate that this may be the point at which the filament W1 feeds the disk.

We posit the presence of a wide angle outflow (half-opening angle of ∼45°) along the same (three-dimensional) axis as the jet. The diffuse and blue filamentary structures on scales of 1.5″ − 4″ to the north participate in this conical radial outflow (and are thus in front of the nucleus). Equivalently, the redshifted gas at the same distance to the southwest is behind the nucleus and part of the redshifted conical outflow. Apart from this velocity structure, the conical outflow model is supported by the velocity gradients seen in the blueshifted (north) and redshifted (south) ionized gas found in the regions 1″ − 4″ from the nucleus. With a radial velocity of ∼400 km s⁻¹, the outflow needs to be powered by the AGN, but since the blue- and redshifted cones are only partially filled with ionized gas, the ionized gas would have to be entrained from the patchy and filamentary interstellar medium by a wide angle wind launched near the black hole. Each cone must be at least 3″ in extent, and thus the crossing time is 0.6 Myr, significantly shorter than the lifetime of the M 87 larger-scale jet and lobes (Owen et al. 2000).

Many of the filaments seen in the inner 4″ appear to cross the nucleus. We see all of the following three cases: (a) those that follow the potential of the black hole, and thus likely truly enter the inner 5 pc; (b) those that enter the nucleus, but only in projection as their bulk blueshifts are not changed as (projected) nuclear separation decreases; and (c) at least one mixed case where a highly blueshifted filament participates in a systemic-velocity-offset quasi-Keplerian rotation as it crosses the nucleus in projection; for illustration, we modeled this as Keplerian rotation in a plane with inclination = 60°. Two potential explanations are (i) a gas filament with a large bulk velocity that approaches the nucleus at a high inclination and (ii) Keplerian rotating gas near the nucleus that is swept into the biconical outflow.

With respect to the inner ionized gas disk: if the M 87 jet is intrinsically two-sided, then the northwestern jet, which is brighter due to Doppler boosting, is oriented toward us. For the inner accretion disk (< 0.1″) to be perpendicular to the jet axis, we expect that the southeastern half of the disk is the near side, and the northwestern the far side. If this is maintained (i.e., if the disk warp is “slight”), we expect the same for the arcsecond-scale ionized gas disk. Indeed, dust lanes seen in the HST images of Ford et al. (1994) led them to suggest that the southeastern side is the near side.

The twisted isophotes seen in the inner arcsecond of the ionized gas disk potentially resolve another long-standing problem: the mismatch between the axes of the inner ionized gas disk and the prominent jet. This issue is extensively discussed in Jeter et al. (2019). Our velocity maps in Fig. 9 provide clear evidence for twisted velocity isophotes in the inner arcsecond of the disk. The twist is clearest from PA ∼ 42° at 0$\stackrel{″}{.}$3–0$\stackrel{″}{.}$5 to PA ∼ 90° at 0$\stackrel{″}{.}$1–0$\stackrel{″}{.}$2. Inside this radius, there is a suggestion of a further twist: the value of this innermost PA is most clearly constrained using the residual maps of our best-fit model (bottom-right panel of Fig. 9). Here, the residual over-subtraction is at PA ∼ 17°, perpendicular to the jet axis. Since our innermost resolved region is about ∼5 pc, there is of course room for further twists down to the jet-launching region around the black hole. For radii ≥0$\stackrel{″}{.}$3, the redshifted side of the disk shows some indications of further twists in the velocity isophotes, though at tens of arcsecond scales we again support a PA close to 40° and an inclination equal to that of the inner arcsecond disk.

Our kinematic analysis clearly shows the presence of multiple velocity components (disk, outflow, and filaments) overlapping in the nucleus. Given this, it is dangerous to constrain the inclination of the disk using isophotal fits. Our Hα+[N II] and [O I] moment 0 maps do not show a clear and consistent inclined disk shape. While the different intensity levels are highly non-axisymmetric, from these maps one gets more of an impression of a face-on, rather than inclined, disk. At least a part of this face-on appearance could be due to filament 8 passing through the nucleus in projection. Its zero velocity gas would extend the intensity contours of the disk along its minor axis. To the north, ionized gas in the outflow is seen even in the inner arcsecond, elongating the north-south extent at some radii. Morphologically, the disk is best constrained via its spiral arm structure, which we see both in intensity and kinematics. These spiral arms are better encompassed by a i = 25° rather than i = 42° disk. The kinematics of the [O I] line in the inner 0$\stackrel{″}{.}$2 to 0$\stackrel{″}{.}$6 is slightly better fit with an inclination of ∼25° as compared to 42°. Finally, we have also shown that the kinematics of the large-scale filaments W1 and W2, which extend from a few to tens of arcseconds from the nucleus, can be reasonably explained in most parts by a rotating disk with i = 25° (but not as well or consistently with an i = 42° disk). Together, the three results favor a relatively face-on disk rather than one at i ∼ 42°. This lower-inclination disk is also consistent with the posited three-dimensional direction of the jet. However, given the warp in the PA of the disk, it is impossible to rule out a change in inclination inside our central resolution elements.

A precise measurement (with well-defined errors) of the black hole mass via ionized gas kinematics is not possible due to the many morphological and kinematic complexities in the ionized gas, which have been discussed extensively in this work. Even removing the effects of the ionized filaments that cross the nucleus (sometimes only in projection), the outflow component, the irregular shape of the disk, the non-axis-symmetric velocities displayed by the disk, the twisted inner velocity isophotes, and the finding that the spiral arms are slightly redshifted from the diffuse gas in the inner disk all preclude a sufficiently accurate determination of the inclination of the subarcsecond-scale disk, and thus the black hole mass.

Nevertheless, the deep and integral spectroscopy from MUSE allows us to clarify and constrain previous contradictory black hole mass measurements: for example, the HBHs of G11, the EHTC, and Liepold et al. (2023) versus the LBH of W13. First and foremost, we have shown that the complexities of the nuclear ionized gas highly complicate the determination of an accurate measurement of the black hole mass from ionized gas kinematics, which immediately lends more weight to the high mass value. However, some indications can be leveraged from the new data, especially via the [O I] velocities. With the caveat of all the complexities involved, we find no reason to favor a i = 42° disk over a lower one (i ∼ 25°). However, several factors support, albeit sometimes weakly, a disk inclination closer to 25°: (a) fits to the subarcsecond ionized gas disk with the Kinemetry package support inclinations closer to 20°–25° and a PA of ∼45°, with the disk velocities dominated by rotation; (b) a minimization of the differences between the observed velocity fields of the subarcsecond ionized gas disk and simulated velocity fields derived via KinMSpy over a range of black hole masses and inclinations lends weak support to the scenario of a HBH in a low-inclination disk, as opposed to a LBH in a higher (42°) inclination disk; (c) visual comparisons of the LBH i42 model and the HBH i25 model with the observed velocity fields and pv diagrams also weakly support the HBH i25 model, though their predictions are very similar within the complexities, and large linewidths, of the emission lines; (d) morphologically, the spiral arms in the inner arcsecond of the ionized disk, which are the most plausible parts of the disk involved in circular rotation, are better encompassed by a i = 25° rather than i = 42° disk; and (e) the outer ionized gas filaments on scales of tens of arcseconds can be relatively well fit with a rotating disk with i = 25°, but not with i = 42°. At these large nuclear distances, the rotation is dominated by the galaxy potential, so there is no longer a degeneracy between the black hole mass and inclination in our rotational fits. Finally, we note that the lower-inclination disk provides a more consistent alignment between the axes of the ionized gas disk, the well-known jet, and the biconical outflow we find here. Thus, the most consistent explanation for the observed inner ionized gas rotating disk is a black hole with a mass closer to 6.0 × 10⁹ M_⊙ surrounded by a low-inclination disk (i ∼ 25°).

To reduce the tension between the G11 and W13 mass measurements, Jeter et al. (2019) suggest intrinsic sub-Keplerian disk rotation of the ionized gas as a potential solution. We tested the specific RIAF model proposed by these authors (over a range of black hole masses and disk inclinations) against the observed velocity field of the subarcsecond nuclear disk and find the following: (a) the sub-Keplerian rotational velocities of RIAF models indeed allow a higher black hole mass in a given disk inclination as compared to Keplerian rotation models, but the specific value of Ω in the Jeter RIAF model is not sufficiently small to allow a HBH in a i = 42° RIAF disk to masquerade as a LBH in a i = 42° Keplerian disk; and (b) the residual (observed minus model) velocity fields of Keplerian models were compared to those with RIAF models. Except for one masking scheme used by us (Gauss-fit masking), Keplerian models result in significantly smaller velocity residuals as compared to RIAF models when both the black hole mass and inclination can freely vary. Reducing the value of Ω can solve the issues in item (a) but not item (b). Furthermore, the disk outflow in the RIAF model should result in significant noncircular motions in the disk. However, our Kinemetry fits to the subarcsecond ionized gas disk imply that the bulk of the velocities can be fit with pure rotation, that is to say, that nonrotative contributions are small. Thus overall, a HBH in a low-inclination disk is a simpler and better-fit alternative to a HBH in a RIAF inflow.

Putting our results together gives us a simple and consistent picture of the nucleus of M 87. A black hole mass closer to ∼6.0−6.5 (rather than 3.5) billion solar masses (for a distance of 16.8 Mpc), together with a warped (at the PA at least) subarcsecond ionized disk with inclination ∼25°. The PA of the innermost accretion disk seen in our residual [O I] velocity map suggests a disk axis aligned with the axes of both the jet and the (partially filled) nuclear radial biconical outflow. The disk inclination to the LOS matches that expected from Doppler boosting of the approaching jet. The kinematics of the ionized filaments at 4″–20″ are reasonably well explained by a (partially filled) rotating disk with inclination 25°, though the several kinks and bubbles along their length point to the presence of other flow velocities and instabilities.

The MUSE WFM and NFM cubes are synergetically useful as they can also be used to extract the stellar dynamics on scales of 0$\stackrel{″}{.}$1 to 30″. Schwarzschild modeling of the observed intensity, velocity, and LOS velocity distribution of the stars will be presented in Osorno et al. (in prep.). Furthermore, the MUSE cubes, together with multi-epoch and multi-filter HST imaging, provide a deep insight into the (patchy) dust structures, the ionized gas distribution, and their relationship with the parsec-scale to kiloparsec-scale radio jet. This is explored in Richtler et al. (in prep.).

Acknowledgments

This work was funded by the National Agency for Research and Development (ANID). J.O. acknowledges support from/Scholarship Program/Doctorado Nacional/2017-21171690. N.N. and J.O. acknowledge funding from ANID via Nucleo Milenio TITANs (NCN19−058), Fondecyt 1221421, and Basal FB210003. P.H. is a member of, and received financial support from, the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. Based on observations collected at the European Southern Observatory under ESO programme 0103.B-0581.

References

Appendix A: VLT/MUSE-NFM data versus HST/STIS data

Previous ionized gas based measurements of the black hole mass in M87 used aperture (Harms et al. 1994) or slit spectra (Macchetto et al. 1997; Walsh et al. 2013) from instruments on HST. In this work, we use MUSE IFU data cubes; the new MUSE data are superior for several reasons. MUSE is a relatively high throughput instrument on a 8 meter class telescope: in collecting area the VLT is a factor of ∼11 larger than HST. In AO mode, MUSE delivers a higher spatial resolution than STIS. The STIS data used in Walsh et al. (2013) consist of spectra in five parallel slits, one of them crossing the nucleus of the galaxy, whereas the MUSE NFM datacube fully samples the inner 8″ × 8″. The larger FOV and full spatial sampling allows a more comprehensive analysis along any PA. The MUSE data also covers a larger wavelength range. The only drawback of the data set is that MUSE has a slightly lower spectral resolution as compared to previous HST instruments used for such studies. However, in a massive elliptical like M87, with relatively large velocities and velocity dispersions, the MUSE spectral resolution is sufficient to accurately resolve and measure gradients in these.

Figure A.1 compares the pv diagrams obtained from the STIS data of Walsh et al. (2013) and our MUSE data along a nuclear slit at PA 51° (the slit PA used by W13). There is a improvement of the quality of the data, especially in the weaker [O I] λ6300 line.

Fig. A.1. Comparison of the pv diagrams obtained from previous HST/STIS long-slit spectra and our VLT/MUSE-NFM-AO datacube for the [N II] λ6583 (top panels) and [O I] λ6300 (bottom panels) emission lines. The pv diagrams are along a slit through the nucleus of M87 at PA = 51°. The horizontal white lines correspond the recessional velocity of M87 adopted by W13: Vrecc = 1335 km s−1.

Appendix B: Software and analysis methods

Section 3 provides a brief overview of the software and methods used in this work. In this appendix, we provide a more detailed description of these.

The stellar kinematics was measured using the version of pPXF (Cappellari & Emsellem 2004; Cappellari 2017) included in the GIST Pipeline (Bittner et al. 2019). The mentioned technique selects a linear combination of stellar templates that best fit the observed galaxy spectrum and delivers stellar moment maps (total intensity, radial velocity, velocity dispersion, and higher moments). Before running pPXF, spaxels were Voronoi binned (Cappellari & Copin 2003), setting a minimal S/N in each bin of 300 for the WFM datacube and 80 for the NFM datacube. The continuum S/N for the Voronoi binning was calculated within the wavelength range 5700 – 5750 Å, and spaxels with S/N of less than 3 in the continuum were excluded from the binning. Bins that included spaxels affected by emission from the jet and its knots were excluded from the analysis. We used the MILES stellar template library (Vazdekis et al. 2010), which covers 3525 – 7500 Å and includes the Mg λ5177 and Na λ5896 lines. The GIST Pipeline was run on the wavelength range 5050 – 6000 Å for the WFM datacube, following Sarzi et al. (2018) and Emsellem et al. (2014); and the wavelength range 5240 – 5770 Å for the NFM datacube, in order to avoid the Mg λ5177 line.

Following the procedure of Sarzi et al. (2018), we chose to use four moments: the velocity (M1), the velocity dispersion (M2), and the third and fourth coefficients of the Gauss-Hermite polynomials, representing the skewness (M3) and kurtosis (M4) of the LOS velocity distribution. We used a sixth (tenth) degree additive polynomial for the NFM (WFM) datacube. We emphasize that GIST Pipeline provides the integrated intensity (M0) per spaxel (not per bin).

The GIST Pipeline includes a module to model the star formation history and derive the weighted mean of the age and metallicity of the templates used to fit the absorption lines. These are also derived via pPXF, but in a different run from that described above for the kinematics. Here, for both WFM and NFM datacubes, we requested that four moments be extracted and used a fourth-degree multiplicative polynomial. The wavelength ranges were the same as used in the stellar kinematics run.

The ionized gas kinematics can be also modeled with GIST Pipeline using GANDALF (Sarzi et al. 2006; Falcón-Barroso et al. 2006; Bittner et al. 2019). The program separates the contribution of the emission lines from the stellar absorption lines in each bin. This step is performed using the results of the pPXF run as inputs; for our analysis, we used the results from Voronoi binning. We chose to fit the wavelength range 4750−6800 Å, since it covers the main emission lines in the MUSE range.

While binning improves the S/N of the spectrum and thus reliability of the resulting moment maps, spatial resolution is lost in key areas of the maps. For this reason, we also used the emission line moment maps that result from a single Gaussian fit to each spectral line in each spaxel’s observed spectrum.

Since nuclear ionized gas velocities in M87 reach almost ∼1000 km s⁻¹, relatively close line complexes, for example Hα+[N II] λλ6548,6583, are blended in the nucleus. We thus analyzed the principle ionized gas lines in groups whose wavelengths are close together. For each group, a wavelength range that covers all the lines plus continuum was selected. A straight line was fit to the continuum in order to obtain an emission-line-only spectrum. The intensity, velocity and dispersion maps of all the lines were obtained from a single Gaussian fit to each line in the group. To facilitate the fitting, two conditions were imposed: (a) all lines in a given group have the same velocity and velocity dispersion; and (b) the amplitude ratios of the following lines were fixed: [O III] λ4959/[O III] λ5007 = 0.35; [O I] λ6364/[O I] λ6300 = 0.333; and [N II] λ6548/[N II] λ6583 = 0.34.

Appendix C: PV diagrams of NFM filaments

Figure C.1 shows the pv diagrams of the NFM filaments presented in Sect. 7.2, but not shown in Fig. 4.

Fig. C.1. Same as Fig. 5 but for the remaining four of the eight filaments marked in the NFM moment maps (dotted black curves in Fig. 4).

Appendix D: Additional masks for the residual velocity maps

Fig. D.1 shows other masks applied to the residual velocity maps and not shown in Fig. 14.

Fig. D.1. Masks applied to the subtracted velocity maps (see also Fig. 14): [O I] λ6300 Gauss-fit mask (top left), [O I] λ6300 spiral-arm mask (top right), [O I] λ6300 annular mask (bottom left), and Hα+[N II] λλ6548,6583 Gauss-fit mask (bottom right).

Appendix E: Biconical outflow models

Figure E.1 shows two models of a biconical outflow, to illustrate alternative scenarios to the biconical outflow discussed in Sect. 11 and shown in Fig. 17. Here we show the velocity field expected if both cones are fully filled with ionized gas. As the observed blueshifted region to the north is wider than the redshifted region to the south, the half opening angle of the receding cone is illustratively chosen to be smaller than the approaching cone: the approaching (receding) cone has a half opening angle of 45° (30°). The model in the left panel creates approximately the observed blue- and redshifted areas, but not the velocity gradients within these, and the outflow axis is perpendicular to the jet axis. In the right panel, the cone axes are chosen to project to the PA of the jet, but the resulting velocity field fails to reproduce the observed velocities.

Fig. E.1. Same as the middle and right panels of Fig. 17, but this time showing both cones superimposed in a single panel. The left panel illustrates a geometry that allows blueshifts to the north and redshifts to the south, as seen in the left panel of Fig. 17. The axes of the cones in this case project to PA 10°, roughly perpendicular to the jet axis. The right panel illustrates the case in which the cone axis is aligned with the jet axis, i.e., at PA −72°.

All Figures

Fig. 1. Maps of stellar properties derived with the GIST Pipeline run on the MUSE WFM datacube. Left: stellar velocity map; zero velocity corresponds to a radial velocity of 1313 km s−1. Middle: map of the weighted mean stellar age of the best-fit templates. Right: map of the weighted mean metallicity of the best-fit templates. These maps were derived using GIST Pipeline over a wavelength range of 5050−6000 Å, with Voronoi binning used to obtain a S/N ≥ 300 in each bin. All spaxels with significant jet emission were masked in the datacube to avoid confusion.

In the text

Fig. 2. Moment and residual maps for the [N II] λ6583 line from the MUSE WFM datacube. Top: moment maps of the [N II] λ6583 line as derived from Gaussian fits to the datacube. From left to right are the total intensity, velocity, and dispersion maps. Two ionized gas filaments (dotted black lines) are indicated by green arrows and the labels W1 and W2; these dotted lines are the pseudo-slits along which the pv diagrams of Fig. 3 are extracted. Several ionized gas bubbles are indicated by purple arrows, and the position of the outflow is indicated by a red arrow. Bottom: residual velocity maps of the [N II] λ6583 line in the datacube. These maps, discussed in Sect. 7.1, show the residual velocity field after subtracting a model of a rotating thin nuclear disk (in a SMBH plus stellar potential) for a black hole and nuclear disk with the following parameters (M•, disk inclination): 6.0 × 109 M⊙ and 42° (left); 3.5 × 109 M⊙ and 42° (middle); and 6.0 × 109 M⊙ and 25° (right; our preferred model; see Sect. 7.1); all have a disk PA of 45°.

In the text

Fig. 3. Position-velocity diagrams of the [N II] λ6583 line for pseudo-slits along the gas filaments W1 and W2 identified in Fig. 2. The x-axis zero position of each filament corresponds to the black cross in Fig. 2. The three curves are the expected disk rotation velocities with the parameters described in Sect. 7.1 and Fig. 2, with colors as listed in the insert. At offsets beyond 10″, the velocities of both filaments are best fit with the i = 25° model. This figure is discussed in Sect. 7.1 but placed here for easy comparison with Fig. 2.

In the text

	Fig. 4. [N II] λ6583 emission line moment maps from Gaussian fits to the NFM cube. From left to right are maps of the total intensity, velocity, and velocity dispersion. Eight ionized gas filaments (dotted black lines) and three pseudo-slits (dotted magenta lines) are numbered and indicated with arrows of the same color in all panels.
In the text

	Fig. 5. Same as Fig. 3 but for four of the eight filaments marked with dotted black curves in the NFM moment maps of Fig. 4. The remaining four filaments are presented in Fig. C.1. These pv diagrams are discussed in Sect. 7.2 but are placed here for easy comparison with Fig. 4.
In the text

	Fig. 6. Same as Fig. 3 but for pseudo-slits C1 and C2, marked with dotted magenta curves in the moment maps of Fig. 4.
In the text

Fig. 7. Extra pv diagrams of circular pseudo-slits in the [N II] λ6548 moment maps (Fig. 4) Left: same as Fig. 3 but for the pseudo-slit C3 (nuclear radius 0$\stackrel{″}{.}$75), which is marked with a dotted magenta line in the left and central panels of Fig. 4. Right: same as the left panel but for a pseudo-slit at a nuclear radius of 0$\stackrel{″}{.}$375. The range of velocities is larger in order to show the full Hα+[N II] line complex.

In the text

	Fig. 8. Zoomed-in view of the NFM moment 0 image of the [N II] λ6583 line to better illustrate the spiral arm pattern and disk morphology. For reference, the blue and red ellipses show the expected morphology of a thin disk with inclinations of 42° and 25°, the values favored by W13 and this work, respectively.
In the text

Fig. 9. Same as Fig. 4 but for the [O I] λ6300 line in the MUSE NFM cube. Only a subset of the NFM FOV, in which the S/N of the [O I] line is high enough, is shown. For illustration, the blue and red ellipses in the top panels show the expected morphology of a thin disk with an inclination of 42° and 25°, the inclination values favored by W13 and this work, respectively. In the rightmost panels, a cross marks the kinematic center; the apertures marked with black circles and labeled with numbers 1 to 6 are used in Fig. 10.

In the text

Fig. 10. [O I] spectral profiles (solid black lines) in selected apertures of the nuclear ionized disk for the apertures indicated in the right panels of Fig. 9, with their single Gaussian fit overlaid in red. The strong line is [O I] λ6300 and the redshifted weaker line is [O I] λ6364. The systemic velocity (dotted black line) and rotation velocities in models with M• = 6.0 × 109 M⊙ and i = 25° (HBH i25; purple) and M• = 3.5 × 109 M⊙ and i = 42° (LBH i42; cyan) are overplotted. Apertures are labeled by their number and their RA and Dec offset from the kinematic center in arcseconds. Apertures were selected to sample portions of the disk with different velocity dispersions and velocity residuals. Aperture 1 is located in the lowest dispersion area of the disk. Aperture 6 shows two components: one centered on the expectation of rotation and the other, with a potentially larger flux, is redshifted 950 km s−1 from the rotating component and 330 km s−1 from the systemic velocity. Apertures to the north (1, 2, and 3) show smaller dispersions, plus a potential component blueshifted by ∼350 km s−1 that is much weaker in flux than the rotating component.

In the text

Fig. 11. Results from Kinemetry for the subarcsecond ionized gas velocity maps. The top two panels show results when both disk PA and inclination are allowed to vary, while the bottom two panels show results when the PA is fixed at 45° and only inclination is allowed to vary. Specifically, the top panels show the variation in the best-fit PA (left) and inclination (right) with nuclear radius. The bottom-left panel shows the variation of the inclination when the PA is fixed to 45°. The velocity maps used were the single Gaussian fit velocity maps of [O I] (red) and Hα+[N II] (green) and the “collapsed” moment 1 map of the [O I] line (black). In each of these three panels, the horizontal lines in the corresponding color show the mean value over radii between between 0$\stackrel{″}{.}$1 and 0$\stackrel{″}{.}$8. The bottom-right panel shows the rotation curve (i.e., the “k1” parameter from Kinemetry; connected blue dots) derived by Kinemetry for the Gaussian-fit [O I] velocity map, and compares this to expectations of Keplerian rotation for the HBH i25 (solid blue line) and LBH i42 (dotted blue line) models, all following the left-hand y-axis. The red line, which follows the right-hand y-axis, shows the k5/k1 ratio (i.e., noncircular velocity to circular velocity ratio) derived by Kinemetry: low values here imply that disk velocities are dominated by rotation. At very small radii, too few pixels are included in each ring and Kinemetry results are unreliable.

In the text

	Fig. 12. Position-velocity diagrams of the Hα+[N II] λλ6548,6583 lines in the WFM cube along the W13 major (45°, left) and minor axes (135°, right). The models (see the inset) and the zero velocity are plotted with respect to the [N II] λ6583 line.
In the text

Fig. 13. Position-velocity diagrams of the Hα+[N II] λλ6548,6583 (left and central) and [O I] (right) emission lines in the NFM cube, along the W13 major (45°, left and right) and minor axes (135°, central). Rotation models overplotted follow the inset. The dotted green line is for an illustrative rotating model with M• = 6.0 × 109 M⊙ (HBH), an inclination of 60°, and a constant blueshift of 650 km s−1 from the systemic velocity.

In the text

	Fig. 14. Examples of masks applied to the subtracted velocity maps. The masks were made by removing pixels with velocity dispersions larger than a given limit (sigma-clip mask). From left to right are the masks for [O I] λ6300 and Hα+[N II] λλ6548,6583. See Fig. D.1 for other types of masks.
In the text

Fig. 15. Standard deviation (in km s−1) of values in the residual (observed minus KinMSpy model) velocity map, as a function of the black hole mass and disk inclination of the model, following the color bar to the right of each panel. The left column is for the Keplerian disk model, and the right column for the Jeter RIAF model. From top to bottom are the standard deviations measured using Gauss-fit, sigma-clip, and spiral-arm masks. The panels of each row follow the same color bar.

In the text

	Fig. 16. Same as Fig. 15 but for the WFM Hα+[N II] λλ6548,6583 velocity map (i.e., the extended filaments) and the Keplerian rotation model. Both panels show the standard deviation of the residual velocity map, with the left panel showing the result with the Gauss-fit mask and the right panel the sigma-clip mask.
In the text

Fig. 17. Modeling of the ionized gas outflow in M 87. Left panel: residual (after subtracting the HBH i25 model for the nuclear disk) [N II] velocity map. We emphasize the gradients in the blue- and redshifted velocities to the northeast and south. Middle and right panels: models of a filled conical outflow made with KinMSpy, with each cone shown separately. The cone half-opening angles are 45° (approaching; middle panel) and 30° (receding; right panel). Each cone has an extent of 3″ along its axis, and the radial outflow velocity within each cone is 400 km s−1. The inclination of both cone axes to the LOS is 18°. The green line shows the projection of the cone axes on the plane of the sky; these project to PA 288°, the PA of the jet axis. The dotted orange lines delineate the projected shape of a disk perpendicular to the cone axes.

In the text

	Fig. A.1. Comparison of the pv diagrams obtained from previous HST/STIS long-slit spectra and our VLT/MUSE-NFM-AO datacube for the [N II] λ6583 (top panels) and [O I] λ6300 (bottom panels) emission lines. The pv diagrams are along a slit through the nucleus of M87 at PA = 51°. The horizontal white lines correspond the recessional velocity of M87 adopted by W13: Vrecc = 1335 km s−1.
In the text

	Fig. C.1. Same as Fig. 5 but for the remaining four of the eight filaments marked in the NFM moment maps (dotted black curves in Fig. 4).
In the text

	Fig. D.1. Masks applied to the subtracted velocity maps (see also Fig. 14): [O I] λ6300 Gauss-fit mask (top left), [O I] λ6300 spiral-arm mask (top right), [O I] λ6300 annular mask (bottom left), and Hα+[N II] λλ6548,6583 Gauss-fit mask (bottom right).
In the text

Fig. E.1. Same as the middle and right panels of Fig. 17, but this time showing both cones superimposed in a single panel. The left panel illustrates a geometry that allows blueshifts to the north and redshifts to the south, as seen in the left panel of Fig. 17. The axes of the cones in this case project to PA 10°, roughly perpendicular to the jet axis. The right panel illustrates the case in which the cone axis is aligned with the jet axis, i.e., at PA −72°.

In the text