11institutetext: Departamento de Astronomía, Universidad de Concepción, Barrio
Universitario S/N, Concepción, Chile 22institutetext: Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany 33institutetext:
Department of Astronomy, University of Texas at Austin, Austin, TX 78712-1205,
USA 44institutetext: Department of Astronomy, University of Michigan, Ann
Arbor, MI 48109, USA

# Revisiting the black hole mass of M87* using VLT/MUSE Adaptive Optics
Integral Field Unit data I.

Ionized gas kinematics

J. Osorno E-mail<EMAIL_ADDRESS>N. Nagar E-mail<EMAIL_ADDRESS>T. Richtler 11 P. Humire 22 K. Gebhardt 33 K. Gultekin 44

(Received March 31, 2023; accepted XXX XX, 2023)

###### Abstract

Context. The stellar dynamic-based black hole mass measurement of M87 is twice
that determined via ionized gas kinematics; these values disagree by more than
3$\sigma$. The former is closer to the mass estimated from the diameter of the
gravitationally-lensed ring around the black hole.

Aims. Using a deeper and more comprehensive ionized gas kinematic dataset, we
aim to better constrain the complex morphology and kinematics of the nuclear
ionized gas, thus gaining insights into the reasons behind the disagreement of
the mass measurements.

Methods. We use new Narrow Field Mode with adaptive optics, and Wide Field
Mode integral field spectroscopic data from the Multi Unit Spectroscopic
Explorer instrument on the Very Large Telescope, to model the morphology and
kinematics of multiple ionized gas emission lines (primarily H$\alpha$+[N  II]
$\lambda\lambda$6548,6583 and [O  I] $\lambda$6300) in the nucleus of M87.

Results. The new deep dataset reveals complexities in the nuclear ionized gas
kinematics which was not seen in earlier sparse and shallower Hubble Space
Telescope spectroscopy. Several ionized gas filaments, some with large flow
velocities, can be traced down into the projected sphere of influence.
However, not all of these truly pass close to the black hole. Additionally, we
find evidence of a partially-filled biconical outflow, aligned with the jet,
with radial velocities up to 400 km s-1. The nuclear rotating ionized gas
’disk’ is well resolved in our datacubes: its twisted velocity isophotes
indicate a warped disk. The complexity of the nuclear morphology and
kinematics (the mix of a warped disk with spiral arms, large linewidths,
strong outflows, and filaments crossing the black hole in projection)
precludes the measurement of an accurate black hole mass from the ionized gas
kinematics. The observed morphology and kinematics of the ionized gas can be
explained with either a 6.6 $\times\ 10^{9}$ $\rm\leavevmode\nobreak\
M_{\odot}$ black hole with a disk inclination of $\sim$25∘, or a 3.5 $\times\
10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$ black hole with a disk inclination
of $\sim$42∘, with more support for the former. The velocity isophotes of the
sub-arcsecond ionized gas disk are twisted and the position angle of the
innermost ($\lesssim$5 pc) gas disk tends towards a value perpendicular to the
radio jet axis.

Conclusions. An inclination close to 25∘ (rather than the previously posited
42∘) for the nuclear gas disk, and the warp in the sub-arcsec ionized gas
disk, help to reconcile the contradictory nature of key earlier results: (a)
the mass discrepancy between stellar and ionized gas black hole masses (our
results support the former); (b) the mis-orientation between the axes of the
ionized gas disk and the jet (we find these aligned in both 2- and
3-dimensions). Further, we identify a previously unknown 400 km s-1
(partially-filled) biconical outflow along the (3-dimensional) jet axis, and
show that the velocities of the two largest ionized gas filaments at 8″–30″
nuclear distances can be explained primarily by rotation in the extension of
the (inclination $\sim$25∘) nuclear ionized gas disk.

###### Key Words.:

galactic dynamics – black hole mass – M87

##  1 Introduction

There is increasing evidence that all the 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 $\sim$200–250 galaxies - are
enabled via molecular (e.g., Barth et al., 2016) and ionized gas kinematics
(e.g., Schnorr Müller et al., 2011), stellar dynamics (e.g., Rusli et al.,
2011) and reverberation mapping (RM; e.g., Bentz et al., 2010). For this
relative 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 RM’ method (using the width of a broad emission line and the
luminosity of the nearby continuum in broad line active galactic nuclei or
AGN) can be used to estimate the SMBH mass. These scaling relationships have
been used to understand the role of the 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 (M87) where 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$\sigma$ of each quoted measurement.

NGC4486 (M87) is a giant elliptical galaxy located in the center of the Virgo
cluster. It has a prominent relativistic jet which has been well studied on
scales of 10 gravitational radii (Event Horizon Telescope Collaboration, 2019)
out to $\sim$40 kpc from the nucleus (Jeter & Broderick, 2021). The jet is
projected on the sky with a position angle (PA) of 288∘ (Walker et al., 2018).
Jeter et al. (2019) argued that the jet has an inclination of 18∘ to the line
of sight and found 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 M87 as a test case for black hole mass
measurement (in)consistency from different methods: it is one of the brightest
galaxies in the sky; as a massive cD elliptical with a large velocity
dispersion, it is expected to host a massive $\geq$
$10^{9}$$\rm\leavevmode\nobreak\ M_{\odot}$ SMBH; and it hosts one of the most
studied and well-known ’radio-loud’ AGN. The first attempts to measure M∙ in
M87 were made by Sargent et al. (1978) and Young et al. (1978), who suggested
a black hole mass of about $10^{9}$$\rm\leavevmode\nobreak\ M_{\odot}$.

In companion papers, Ford et al. (1994) presented H$\alpha$+[N  II]
$\lambda\lambda$6548,6583 emission-line images of the M87 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 M87 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 in PA $\sim$1∘–9∘,
inclination (i) $\sim$42∘, and a black hole mass of 2.4 $\times$ 109
$\rm\leavevmode\nobreak\ M_{\odot}$ assuming a distance of 15 Mpc to M87.
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
$\times$109 $\rm\leavevmode\nobreak\ M_{\odot}$ 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 modelling of observed stellar brightness and stellar
line of sight velocities to measure a black hole mass of (6.6 $\pm$ 0.4)
$\times\ 10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$ for a distance of 17.9
Mpc. They used data from the Integral Field Spectrograph (NIFS) on the Gemini
North Telescope with adaptive optics 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 M87 via the kinematics of
the H$\alpha$+[N  II] $\lambda\lambda$6548,6583 emission lines using data from
the Space Telescope Imaging Spectrograph (STIS) aboard HST. Their spectral
data were taken in five parallel slits in PA = 51∘, separated by
0$\aas@@fstack{\prime\prime}$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$\alpha$ and [N  II] lines are highly blended, they derived a M∙ of
3.5${}^{+0.9}_{-0.7}$ $\times\ 10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$ for
a galaxy distance of 17.9 Mpc.

Recently, the black hole mass of M87 was constrained using data from the 2017
campaign of the Event Horizon Telescope Collaboration (EHT Collaboration,
2019). A SMBH enveloped in radio emitting plasma is expected to be have a
bright gravitationally lensed ring of emission surrounding the black hole
’shadow’. The diameter of ring is predicted to be $\sim$9.6–10.4 GM∙/$c^{2}$,
relatively independent of the SMBH spin (Johannsen & Psaltis, 2010; Gralla et
al., 2020). The ring diameter measured by the EHT, together with standard
General Relativity, imply a M∙ of 6.5 $\pm$ 0.2 (statistic) $\pm$ 0.7
(systematic) $\times\ 10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$ for a
distance to the galaxy of 16.8${}^{+0.75}_{-0.66}$ Mpc. Scaled to this
distance, the G11 and W13 SMBH masses are, respectively,
6.14${}^{+1.07}_{-0.62}$ $\times\ 10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$
and 3.45${}^{+0.85}_{-0.26}$ $\times\ 10^{9}$ $\rm\leavevmode\nobreak\
M_{\odot}$ (Event Horizon Telescope Collaboration, 2019). The EHT result is
thus compatible with G11 but not with W13.

Liepold et al. (2023) studied the stellar kinematics of M87 with data from the
integral field spectroscopic 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 $\sim$5 kpc and outwards).
They derived a SBHM of 5.37${}^{+0.37}_{-0.25}\times
10^{9}$$\rm\leavevmode\nobreak\ M_{\odot}$ (for a galaxy distance of 16.8
Mpc), using Schwarzschild orbit modelling with a varying mass-to-light ratio
(M/L). This result is most consistent with the EHT value, and also relatively
consistent with the G11 value.

The morphology of the ionized gas in M87, on arcseconds to tens of arcseconds
scales, 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 $\sim 3$ kpc to the NW and $\sim 8$ kpc to
the SE. 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 arcmin, and showed
that the filaments have perturbed kinematics with velocity differences of
$\sim$700 – 800 km s-1, but a fairly uniform velocity dispersion of $\sim$100
km s-1 over the filaments.

As mentioned above, the EHT Collaboration adopted a distance to M87 of 16.8
Mpc (thus a scale of 0.081 kpc per arcsec), while G11 and W13 used a distance
of 17.9 Mpc (a scale of 0.087 kpc per arcsec). The heliocentric recessional
velocity adopted by NED for M87 is 1284 km s-1 (from stars). The radial
velocity listed in SIMBAD is 1256 km s-1. The best fit of W13 to the ionized
gas kinematics in HST/STIS data yielded a radial velocity of 1335 km s-1. 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-1. Our equivalent fit
to the full field of view (FOV) of the Narrow Field Mode (NFM) cube yields a
radial velocity of 1301 km s-1. 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-1. In this work, unless explicitly mentioned otherwise, we use
’zero’ radial velocities of 1313 km s-1 and 1260 km s-1 for the stellar and
ionized gas kinematics, respectively. For easiest comparison with G11 and W13,
we use a distance of 17.9 Mpc.

To better understand the difference between the M87 SMBH masses measured using
ionized gas, stellar dynamics, and EHT ring-size, we have obtained new
observations of M87 with MUSE in its adaptive optics (AO) 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 better signal-to-noise
ratio (SNR) than all previous HST/STIS data. Further, unlike the data from
STIS, the MUSE NFM cube fully samples the nuclear 8$\arcsec\times$8$\arcsec$
region and provides a larger wavelength coverage (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 M87
over the central 1′, allow constraints on the larger scale morphology and
kinematics/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 velocity line of sight 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.).

The paper is distributed in the following sections: observations and data
processing in Sect. 2, with emphasis on the NFM cube; our programs and
modelling procedure in Sect. 3; relevant stellar kinematics and properties in
Sect. 4; the ionized gas morphology and kinematics, together with our results
on constraining the black hole mass and outflow properties in Sect. 5; and
finally, the discussion and conclusions of this study in Sect. 6.

##  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 $\AA$ with a wavelength
sampling of 1.25 $\AA$ per pixel in both WFM and NFM modes. The spectral
resolution is R = 3026 at 7000Å. The WFM mode has a FOV of
60$\arcsec\times$60$\arcsec$, and a spatial sampling of
0$\aas@@fstack{\prime\prime}$2 per pixel. The NFM mode has a FOV of
8$\arcsec\times$8$\arcsec$ with a spatial sampling of
0$\aas@@fstack{\prime\prime}$0252 per pixel. Since NFM data are obtained in
adaptive optics (AO) mode with Sodium laser artificial guide stars, no data
are obtained in the wavelength range of 5780 – 6050 $\AA$.

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 DIMM seeing of
0$\aas@@fstack{\prime\prime}$9\. The NFM mode observations were made on
February 2, 2020, as part of project ID 0103.B-0581 (P.I. Nagar). Data were
taken at airmasses between 1.25 and 1.5, with DIMM seeing values of
0$\aas@@fstack{\prime\prime}$4 to 0$\aas@@fstack{\prime\prime}$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 seconds, 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
seconds of on-sky exposure, were downloaded from this archive. A fourth
observation, which was aborted after 83 seconds 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 N) and combined into a single 6300 second 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 M87 is relatively
northern for the VLT and has to be observed at relatively large airmass, the
FWHM of the PSF core is expected to be better than 100 mas (see the MUSE
Exposure Time
Calculator111https://www.eso.org/observing/etc/bin/gen/form?INS.NAME=MUSE+INS.MODE=swspectr
and MUSE User Manual
222https://www.eso.org/sci/facilities/paranal/instruments/muse/doc.html). With
no stars in the FOV, we were unable to measure the true PSF in our dataset.
However, three globular clusters in M87 can be distinguished in a continuum
image of the NFM cube after subtracting a model galaxy bulge. The FWHM 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$\arcsec\times$8$\arcsec$, 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. Appendix A compares the two
datasets, with Fig. 14 highlighting the large improvement the new dataset
brings.

##  3 Methods and software

Figure 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 $\AA$, with Voronoi binning used to obtain a
SNR $\geq$ 300 in each bin. All spaxels with significant jet emission were
masked in the datacube to avoid confusion.

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.

All codes used, both publicly available and our own, are in Python, and we
extensively use 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 signal to
noise ratio of the spectra, and spaxels which include the jet were not
included in the analysis. For details, see 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 results, were used to check our results
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, e.g.,
H$\alpha$+[N  II] $\lambda\lambda$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 build position-velocity (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 Sect. 5.3, we assume that a component of the ionized gas is rotating in a
thin disk, and model the rotation, 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 derived using:

$v_{gas}=\sqrt{\frac{G\left[M_{\bullet}+M_{bulge}(r)\right]}{r}}$ (1)

Where $r$ is the nuclear distance, $G$ the gravitational constant,
$M_{\bullet}$ 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:

$M_{bulge}(r)=\int_{0}^{r}4\pi r^{\prime
2}\Upsilon(r^{\prime})\nu(r^{\prime})dr^{\prime}$ (2)

Where $\Upsilon$ is the radially varying mass to light ratio (M/L) of the
stars, $\nu$ is the volume density luminosity of the stars, and $r^{\prime}$
is the integration variable. The V-band luminosity density of M87, as a
function of radius, is plotted in Fig. 1 of Gebhardt & Thomas (2009). The
values in this plot were extracted with version 3.5.1 of the Plot Digitizer
Online Application333https://plotdigitizer.com. 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 $\gtrsim$ 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 e.g.
values of V-band M/L of 8.65 and 3.46 in the central and the outer part of
M87, respectively.

The line-of-sight (LOS) velocity of the ionized gas in a given spaxel or along
a slit was obtained by projecting $v_{gas}$ with the following equation:

$v_{los}=v_{gas}\sin{i}\cos{\theta}$ (3)

Where $i$ is the inclination of the disk respect to the plane of the sky, and
$\theta$ is the angle (in the plane of the disk) between the major axis of the
disk and the deprojected position angle of the slit.

In Sect. 5.4 we discuss the nuclear ionized gas outflow. This outflow is
modelled with one or two cones with vertices in the nucleus. We define our
coordinate system with the $x$ axis positive to the E, the $y$ axis positive
to the N, the observer along the positive z axis, and the origin at the
nucleus. The opening angle, length of each cone, and the radial outflow
velocity, are free parameters. This outflow model is used as input to the
KinMSpy software (Davis et al., 2013), using the feature ’inclouds’. KinMSpy
creates a synthetic data cube that simulates an observation of the outflow
with the characteristics of the observed datacube, which is in turn used to
obtain the simulated velocity maps.

##  4 Stellar Properties from the WFM datacube

A full analysis of the stellar dynamics and Schwarzschild modelling of the
line of sight 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-1 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 $\pm$5 km
s-1 in the inner $\sim$18″with redshift to the NE, and the other of $\pm$20 km
s-1 (in the opposite direction, i.e. redshift to the SW) and going out to
beyond the FOV of the cube. The velocity dispersion map (not shown) is
symmetric with values $\sim$250 km s-1 in a ring of nuclear radius $\sim$20″,
and steadily increasing to more than 400 km s-1 with decreasing nuclear
distance. Given that, as expected, stellar dispersion dominates stellar
rotation in the nucleus; therefore, we will ignore these $\pm$5 or $\pm$20 km
s-1 of stellar rotational velocities when modelling the gas kinematics in the
next sections.

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 $\sim$12″, with values between 13–13.5 Gyr and
0.32–0.35 respectively (except in the innermost arcsec, where the large
dispersion, and the nuclear jet emission, potentially confuse the fits).
Outside $\sim$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 $\sim$23″ (values between 0.2–0.24) and then
slightly increases to values between 0.24–0.3. Note the steep change in the
metallicity at radius $\sim$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$\aas@@fstack{\prime\prime}$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 arcsec of the ionized gas disk, and is mainly relevant in
the analysis of the kinematics of ionized filaments at $\gtrsim$ 2 ″.

##  5 Ionized gas kinematics

Figure 2: Top: moment maps for the [N  II] $\lambda$6583 line as derived from
Gaussian fits to the MUSE WFM datacube. Left to right are the total intensity,
velocity, and dispersion maps. Two ionized gas filaments (black dotted 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] $\lambda$6583 line in the MUSE WFM datacube. These maps, discussed in
Sect. 5.2.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 position angle, disk
inclination): left: 6.6 $\times\ 10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$
(G11), 45∘ and 42∘; middle: 3.5 $\times\ 10^{9}$ $\rm\leavevmode\nobreak\
M_{\odot}$ (W13), 45∘ and 42∘; and right: 6.6 $\times\ 10^{9}$
$\rm\leavevmode\nobreak\ M_{\odot}$ (G11), 40∘ and 25∘ (see Sect. 5.2.1).

Figure 3: Position-velocity diagrams of the [N  II] $\lambda$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. 5.2.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. 5.2.1, but placed here for easy
comparison with Fig. 2.

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. Apart from the rotating disk (with kinematically distinguishable
spiral arms), we see evidence of a nuclear biconical outflow and multiple
filaments at different systemic velocities, some of which cross the nucleus
only in projection. The projection of these multiple components against the
nuclear disk preclude an immediate analysis of the rotating disk. We thus
extensively discuss the ionized gas morphology and kinematics in the following
sub-sections, before discussing the ionized gas disk.

We use both moment maps and pv diagrams of several emission lines covered in
the wavelength range of MUSE, to identify and separate the ionized gas
components in the galaxy. Since the [N  II] $\lambda$6583 emission line is the
brightest in the nucleus, we primarily base our analysis on this line.
However, in the innermost $\sim 1\aas@@fstack{\prime\prime}$5 arcsec, the
lines are very broad and have multiple velocity components, so we are unable
to clearly separate the H$\alpha$ and [N  II] lines in our Gaussian fits. We
thus also rely on the [O  I] $\lambda$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 arcsec.

###  5.1 Moment maps

Ionized gas maps of the central arcmin of M87 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. Figures 2 (top row) shows the
intensity, velocity, and dispersion maps of [N  II] $\lambda$6583 obtained
from the WFM data cubes. The velocity maps assume a systemic velocity of
$1260$ km s-1, and a sigma clipping was applied to remove bad pixels. Several
features are worth noting in these moment maps. The two black dotted lines,
indicated with green arrows, trace the two primary ionized gas filaments which
enter the nucleus (in projection). The filament labelled ”W1” is the most
prominent nuclear filament and extends beyond the MUSE WFM FOV to the W
(Gavazzi et al., 2000; Werner et al., 2010; Boselli et al., 2019). The
filament marked ”W2” also crosses the nucleus but from the E, and connects to
filament W1 about 10″ to the NE 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 $\sim$6″ diameter region to the N of the nucleus
(indicated with the red arrow in the velocity and dispersion maps) that is
$\sim$200 km s-1 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 $\sim$400 km s-1 to the NW (SE).

The dispersion map shows high dispersion to the NE 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
SW 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, i.e. the filament kinematics
does not follow rotation in the potential of the black hole.

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

Figure 5: As in Fig. 3 but for four of the eight filaments marked with black
dotted curves in the NFM moment maps of Fig. 4. The remaining four filaments
are presented in Fig. 15. These pv diagrams are discussed in Sect. 5.2.2, but
are placed here for easy comparison with Fig. 4.

Figure 6: As in Fig. 3 but for pseudo-slits C1 and C2 marked in magenta dotted
curves in the moment maps of Fig. 4.

Figure 7: Left: as in Fig. 3 but for the pseudo-slit C3 (nuclear radius
0$\aas@@fstack{\prime\prime}$75″) marked with a magenta dotted line in the
left and central panels of Fig. 4. Right: as in the left panel, but for a
circular pseudo-slit at a nuclear radius of 0$\aas@@fstack{\prime\prime}$375;
a larger range of velocities is shown in order to show the full H$\alpha$+[N
II] line complex.

The equivalent NFM moment maps are presented in Fig. 4, with the locations of
several pseudo-slits identified. The black dotted 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 magenta dotted
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
N, 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 will be more extensively discussed in Sects. 5.2.2 and
6. At this point we note the diffuse and constant-velocity blue- (red-)
shifted emission to the N (S) of the posited rotating accretion disk which
have high velocities but low dispersion, the multiple filaments to the N (with
velocities apparently lower than the diffuse blue emission), and the high
dispersion in the posited accretion disk. The latter is caused, in part, from
the mixing of multiple velocity components: the disk and outflow (N), and the
disk and a filament crossing the nucleus in projection (S).

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$\alpha$+[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] $\lambda$6300 line as it is relatively
isolated and has sufficient SNR. Figure 8 (top row) shows the moment maps for
the [O  I] $\lambda$6300 line in the NFM data cube; the FOV is smaller here as
we only show regions with relatively high SNR. 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), respectively.

As in the case of the [N  II] line, the primary velocity pattern traces
rotation in a disk, though there are complexities. A clear velocity isophote
twist is seen in the inner disk, the morphology of the outer disk does not
trace a perfect disk shape (and thus cannot be used to precisely constrain the
disk inclination), and the amplitude of the velocity field is not symmetric
across the minor axis. The blueshifted emission has smaller amplitudes than
the redshifted emission in the radius range
0$\aas@@fstack{\prime\prime}$1–0$\aas@@fstack{\prime\prime}$3\. In the next
subsection we show that this is due to a filament which 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
dispersion is also asymmetric: the inner disk has large dispersions, reaching
values of $\sim$600 km s-1, to the SW, and the outer disk has larger
dispersions to the W.

Figure 8: As in Fig. 4 but for the [O  I] $\lambda$6300 line in the MUSE NFM
cube. Only a subset of the NFM FOV, in which the SNR of the [O  I] line is
high enough, is shown. For illustration, the blue (red) ellipses in the top
panels show the expected morphology of a thin disk with inclination 42∘ (25∘),
the inclination values favored by W13 and this work, respectively.  Figure 9:
As in the top left panel of Fig. 8, but for the [N  II] $\lambda$6583 line.
The ellipses here have a smaller semi-major axis to allow a better comparison
to the spiral arms.

###  5.2 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.

####  5.2.1 Wide Field Mode filaments

In Sect. 5.2.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 with the following
parameters:

  * •

Magenta line (hereafter 3.5i42): M∙ = 3.5 $\times\ 10^{9}$
$\rm\leavevmode\nobreak\ M_{\odot}$ (W13), i = 42∘and PA = 45∘.

  * •

White line (hereafter 6.6i42): M∙ = 6.6 $\times\ 10^{9}$
$\rm\leavevmode\nobreak\ M_{\odot}$ (G11), i = 42∘and PA = 45∘.

  * •

Green line (hereafter 6.6i25): M∙ = 6.6 $\times\ 10^{9}$
$\rm\leavevmode\nobreak\ M_{\odot}$ (G11), i = 25∘and PA = 40∘.

The galaxy potential is obtained as described in Sect. 3, and becomes
significant only at radii beyond $\sim$2″. The 3.5i42 and 6.6i42 (magenta and
white) models use the inclination and position angle found by W13 for the
ionized gas disk. The green line is derived using the SMBH mass from G11, but
with values of inclination and position angle modified to better fit our data
(see Sect. 5.3).

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 in 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 $\pm$100 km s-1 (with an extreme value of $\sim$400 km
s-1) over the first $\sim$30″ of the filament. Between 20″ and 30″ there are
two - potentially bubble-shaped - clouds at 90 km s-1 and $-$350 km s-1. 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 line
of sight) or caused by a kink or instability in the filament. From 36″ onward,
the observed velocities follow the overall pattern expected from the 6.6i25
rotating disk model in both values and shape. To make the 6.6i42 and 3.5i42
models agree with the observed velocities over this range, one would require
to change the ionized gas ’zero’ velocity by $\sim$80 km s-1 to the blue,
which would imply an even greater difference between the recessional
velocities of ionized gas and stars ($\sim$1300 km s-1).

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 6.6i25 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 E of the nucleus (Fig. 2),
though we find no explanation for this velocity pattern.

Except in the innermost $\sim$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 6.6i42 (white lines) and 3.5i42 (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 arcsec of the ionized gas disk no
longer plagues us.

####  5.2.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]
$\lambda$6583 line along these filaments and pseudo-slits are shown in Figs.
5, 6, and 15. 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 15) 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 W of the highly blueshifted region
(see 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. 15)
shows differences of up to $\sim$400 km s-1 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 NE of the nucleus. They are the most blueshifted filaments in the inner
nucleus, with blueshifts of $\sim$250–350km s-1 from systemic, and more than
500 km s-1 offset from the rotation model predictions. Nevertheless, their
dispersion velocities are $\leq$ 150 km s-1 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-1 and $-$150
km s-1, they are less blue shifted than filaments 4 and 5. Filament 6 shows a
higher dispersion ($\sim$200 km s-1) 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.1$\arcsec$–0.2$\arcsec$ and
$\sim$1.8$\arcsec$–2.2$\arcsec$. Its dispersion is also small ($\leq$150 km
s-1). 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 $\sim$200 km s-1 from systemic
except over offsets of $\sim$1$\arcsec$–3$\arcsec$ along the pseudo-slit
(effectively the NE region) where a constant blueshift of $\sim$300 km s-1 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 disc, travels to
the SE along a small filament (see 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 arcsec, and
the fact that the kinematics roughly follow the 3.5i42 and 6.6i25 rotating
disk models in the inner two arcsec. We see a small velocity discontinuity at
$\sim$2.5$\arcsec$–3$\arcsec$ offset along the pseudo-slit, which delineates
the separation point between the rotating (inner) and non rotating (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$\aas@@fstack{\prime\prime}$75
centered on the nucleus. The $x$ axis starts (offset = 0) W of the nucleus and
moves clockwise. Over all of the pseudo-slit except the NW quadrant the
observed velocities are well fit with either the 3.3i42 model, or a 6.6i25
model. In the NW quadrant, the observed velocities are systematically
$\sim$200 km s-1 blueshifted from the model predictions. While filaments 6 and
7 of the NFM can be seen as separate components at even larger blueshifts, the
$\sim$200 km s-1 blueshift in the NW does not appear to be an additional
component but rather a warp or systematic bulk velocity in the disk. At this
radius, the S, especially the SW, 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$\aas@@fstack{\prime\prime}$375). Here, the disk velocities
are more symmetric in all PAs, with the 3.3i42 and the 6.6i25 both providing a
reasonable fit. Note that the N blueshifted filaments are still clearly
distinguishable. At this smaller radius, the gas to the SW shows a second
velocity component: a non-rotating component $\sim$300 km s-1 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$\aas@@fstack{\prime\prime}$1–0$\aas@@fstack{\prime\prime}$3\. The
velocity of this non-rotating component matches that of filaments 6 and 7 as
they approach the nucleus in projection.

###  5.3 Ionized gas rotation: a warped low inclination disk

Figure 10: PV diagrams of the H$\alpha$+[N  II] $\lambda\lambda$6548,6583
lines in the WFM cube along the W13 major (45∘, left) and minor axes (135∘,
right). The models (see inset) and the ’zero’ velocity are plotted with
respect to the [N  II] $\lambda$6583 line.

Figure 11: PV diagrams of the H$\alpha$+[N  II] $\lambda\lambda$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). Rotating
models overplotted follow the inset. The dotted green line is for a rotating
model with M∙ = 6.6 $\times$ 109 $\rm\leavevmode\nobreak\ M_{\odot}$, an
inclination of 60∘, and blueshifted by 650 km s-1.

Figure 10 shows the pv diagrams of the H$\alpha$+[N  II]
$\lambda\lambda$6548,6583 line along nuclear slits in 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. Figure 11
shows the equivalent pv diagrams in the NFM cube (left and middle panels)
together with the pv diagram of the [O  I] $\lambda$6300 line along a nuclear
slit in PA = 45∘. In the left panel, the gas
$\sim$1$\aas@@fstack{\prime\prime}$5–4″to the NE of the nucleus (positive
offsets in the figure) shows 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$\aas@@fstack{\prime\prime}$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$\aas@@fstack{\prime\prime}$8 in projection) and then further increases until
it is indistinguishable from gas in the rotating disk. Adding a simple bulk
blueshift to our rotation models cannot reproduce this velocity pattern (the
observed velocities drop too rapidly with increassing 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 3
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-1. 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. 5.4).

Comparing the observed pv data to the models, the 6.6i42 overpredicts the
observed velocities, while both the 3.5i42 and the 6.6i25 models are
relatively indistinguishable when compared to the data. Along the minor axis
(central panel of Fig. 11), the blueshifted gas to the NW is clearly visible
at an almost constant $\sim$200 km s-1 blueshifted velocity. The [O  I] pv
diagram (right panel of Fig. 11) covers a smaller range of offsets. Here the
6.6i42 model is clearly in the upper envelope of the observed velocities.

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 use the [O  I] $\lambda$6300 emission line,
which, while weaker, is relatively isolated from other emission and absorption
lines, and is detected in over the central two arcseconds.

We thus focus in the [O  I] $\lambda$6300 moment maps shown in Fig. 8. The
intensity (moment 0) map shows decreasing flux with increasing nuclear
distance but does not clearly delineate an axisymmetric inclined disk. The
velocity map shows a clear rotational pattern, with a twist close to the
nucleus. In the highest SNR inner ($\leq$0$\aas@@fstack{\prime\prime}$2) disk,
the blueshifted side of the disk has a larger dispersion. Spectra here are
double-peaked with a red component at velocity similar to the redshifted side
of the disk, and a blue rotating component. A high dispersion region
0$\aas@@fstack{\prime\prime}$5 to the W of the nucleus, in an area with a
lower SNR, is also seen: here the extra width could be attributable to one of
the filaments discussed in the previous section.

The twisted ionized gas velocity isophotes in the innermost
0$\aas@@fstack{\prime\prime}$1 (more clearly discerned in Fig. 8, but also
clear on close inspection of Fig. 4) is an important result of this work. At a
radius of 0$\aas@@fstack{\prime\prime}$3, the PA of the disk appears to be
$\sim$45∘ as found by W13. By $\sim$0$\aas@@fstack{\prime\prime}$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. 8). Here the residual map implies a rotation major
axis of $\sim$18∘, which is perpendicular to the PA of the inner jet of M87,
as would be expected if the jet was launched perpendicular to the accretion
disk. Note that Macchetto et al. (1997) also noted the velocity field twist
inside 0$\aas@@fstack{\prime\prime}$2 when modelling their HST spectra with
thin disk models. The twist appears to continue outside the radius of
0$\aas@@fstack{\prime\prime}$3\. Both the [N  II] and [O  I] NFM velocity maps
appear to support a major axis PA closer to $\sim$20$\arcsec$–25$\arcsec$ at a
radius of 0$\aas@@fstack{\prime\prime}$5\. However, the filaments velocity
field in the WFM maps at r $\sim$10s of arcseconds can still be fit well with
a rotating disk in PA $\sim$42∘.

Emsellem et al. (2014) have shown that the stellar velocity field in the
central $\sim$15″ of M87 is twisted, with the stellar kinematic PA changing
rapidly from 17∘ to $-$124∘. In the innermost $\sim$5″, the stars rotate with
a major axis close to PA = 20∘, with the NE side redshifted. This pattern can
also be seen in our WFM stellar map (Fig. 2), where visually the inner
$\sim$10″ stellar kinematics support a PA of 30∘–45∘, with the NE receding,
i.e., similar to the pattern seen in the ionized gas here. The main difference
is that the stellar rotation velocity amplitudes are relatively low
($\lesssim$20 km s-1) and that the twisted isophotes are seen on much larger
scales.

To directly compare the contradictory black hole mass values of G11 and W13,
we construct residual velocity maps of the [O  I] line, by subtracting our
3.5i42 and 6.6i42 rotating disk models (Sect. 5.2.1). These residual velocity
maps are shown in the left (G11) and central (W13) panels in the bottom row of
Fig. 8. For these fixed values of PA and inclination (determined by W13) the
6.6 $\times\ 10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$ black hole
significantly over predicts the LOS observed velocity field, while the 3.5
$\times\ 10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$ black hole shows
relatively smaller velocity residuals, even though it still over predicts the
observed velocity field.

To better explore the parameter space of PA and inclination, we fixed the
black hole mass to 6.6 $\times\ 10^{9}$ $\rm\leavevmode\nobreak\ M_{\odot}$
and varied the PA and inclination to obtain, visually and statistically, a
better fit to the observed velocity field. The best fit is obtained for a PA
of 40∘ and inclination of 25∘. The velocity residual for a model with these
parameters, the 6.6i25 model, is shown in the right panel of Fig. 8. This
residual shows smaller residual velocities for the diffuse ionized gas in the
disk as compared to the left and middle panels. The main residual pattern is a
three or four armed spiral pattern at $\sim$100 km s-1 redshift, embedded in a
blue (outflow or filaments; see the previous subsection) halo. These velocity-
traced spiral arms closely match the spiral arms noted by Ford et al. (1994)
in their HST ionized gas maps, and seen clearly in our moment 0 maps (e.g.,
Fig. 9). The spiral pattern is reminiscent of the ionized gas spirals around
SgrA*, which extend out to $\sim$3 pc, and appear to trace Keplerian circular
or elliptical orbits around SgrA* (Zhao et al., 2009). The innermost
0$\aas@@fstack{\prime\prime}$1 (9 pc) shows an oversubtraction of the model in
PA $\sim$18∘. Rather than a mismatch of the parameters of black hole mass, PA
and inclination, this is most likely the effect of the twisted velocity field
as pointed out in the previous paragraphs. A comparison of the inner disk
residual velocities seen in the middle and right panel of Fig. 8 is shown in
the left panel of Fig. 12. To avoid low SNR regions and the central strong
warp, we use only disk semi-major axes between 0$\aas@@fstack{\prime\prime}$2
and 0$\aas@@fstack{\prime\prime}$5 with a disk inclination of 42∘ or 25∘.
Further, a part of the SW disk, where the derived velocities are confused by
the non-rotating filamentary gas, is masked (see the right panel of Fig. 12).
The 6.6i25 model residuals are more symmetrical and have a smaller dispersion
(47 vs 56 km s-1) compared to the 3.5i42 model residuals, when the masking
uses a 42∘ inclination disk (Fig. 12). In the case of masking with a 25∘
inclination disk, the difference is 53 vs 63 km s-1, respectively. When
varying the inner and outer radius of the mask, the 6.6i25 model residuals are
almost always more symmetrical and with smaller dispersion.

Given this (even if slight) result supporting a lower inclination disk in the
inner arcsecond, we revisit the pv diagrams of Figs. 3, 5, 6 and 7. In these,
the green line corresponds to the G11 value of black hole mass, but the PA and
inclination are now those of our best fit above: 40∘ and 25∘ respectively. The
predicted velocities from this are not very different from the prediction of
the W13 model (magenta line). While there is a partial degeneracy between the
black hole mass and the inclination at these small nuclear distances, the
better fit of the 6.6i25 model to the [O  I] velocity field makes it the more
favored one.

As a further test, we revisit (visually) the morphological inclination of the
ionized gas disk. Ford et al. (1994) used isophotal fits to the inner arcsec
ionized disk (which we know to be a mix of different components, not only a
disk) to derive a PA of 1∘–9∘ (significantly different from the kinematic PA
of $\sim$40–42∘ supported by all kinematic studies which allowed this PA to
vary) and an inclination of 42∘. The [O  I] moment maps (Fig. 8) provide
useful constraints on the disk kinematics, but the moment 0 map is not optimal
for deriving a morphological inclination. We thus show an equivalent map of
the [N  II] $\lambda$6583 line in Fig. 9. Here, 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. 9, the outer ends of the spiral arms are better enclosed in an i = 25∘,
rather than i = 42∘, disk.

We also test if the larger scale filaments seen in the WFM can be explained by
rotation in a (partially filled) ionized gas disk. Here it is important to
note that we use a recession velocity of 1260 km s-1, which is the velocity
which best fits the inner ionized gas disk, and which is 75 km s-1 bluer than
the value used by W13. With our recessional velocity of 1260 km s-1, 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-1
(redshifted). The bottom panels of Fig. 2 show the residual velocity maps of
the [N  II] $\lambda$6583 line after subtracting the 3.5i42, 6.6i42, and
6.6i25 models. At the relatively large nuclear distances seen in the WFM, the
galaxy potential dominates the black hole potential. Both the 6.6i42 and
3.5i42 models (i.e. i = 42∘ models) oversubtract the velocities over most of
the filaments W1 and W2. The i = 25∘ model however produces a velocity
residual closer to zero. As discussed in Sect. 5.2.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 two filaments, both appear to be primarily gas in a rotating disk.

Figure 12: Left panel: histograms of the residual velocity field of the [O  I]
$\lambda$6364 line when using the 3.5i42 (blue) and 6.6i25 (orange) models.
The histograms were calculated over a region as seen in the right panel:
between semi-major radii of 0$\aas@@fstack{\prime\prime}$2 and
0$\aas@@fstack{\prime\prime}$4 for a disk with inclination 42∘, and masking
the high dispersion region to the SW. The residuals of the 6.6i25 model is
more Gaussian-like and has a smaller dispersion (47 vs. 56 km s-1) as compared
to the residuals from the 3.5i42 model.

###  5.4 Ionized gas outflow

The moment maps and pv diagrams in the diffuse blue region and filaments to
the N of the nucleus clearly demonstrate the presence of a non-rotating
component. Additionally, some excessively redshifted gas is seen $\sim$3″ SW
of the nucleus. With projected velocities up to 400 km s-1, the most viable
explanation for these velocities is an AGN driven outflow.

The left panel of Fig. 13 shows the residual velocity map (observed velocity
minus our 6.6i25 rotating disk model) of the [N  II] $\lambda$6583 line, and
emphasizes better not only the blue and redshifted regions to the N and S,
respectively, but also the velocity gradients seen in these.

One could attempt to reproduce these velocity features with a (fully-filled)
biconical outflow. Illustrative examples of these are shown in Fig. 16, 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.
13 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. 13). 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 N (outside
the accretion disk and within the NFM FOV) is primarily between the nucleus
and the observer, and the ionized gas to the SW 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. 13.

Figure 13: Left panel: The residual (after subtracting the 6.6i25 model for
the nuclear disk) [N  II] velocity map. Note the gradients in the blue and
redshifted velocities to the NE and S. Middle and right panels: models of a
filled conical outflow made with KinMSpy, with each cone showed 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 line of sight 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 orange dotted lines delineate the projected shape
of a disk perpendicular the cone axes.

##  6 Discussion and Conclusions

Our new deep MUSE-AO NFM cube, together with the archival MUSE WFM cube, allow
improved constraints on the morphology and kinematics of the ionized gas in
the inner 0$\aas@@fstack{\prime\prime}$1 to 30 arcsec in M87. On scales of
1.5$\arcsec$–4$\arcsec$, we posit that the diffuse and blue filamentary
structures to the N participate in a conical radial outflow (and are thus in
front of the nucleus). Equivalently, the redshifted gas at the same distance
to the SW is behind the nucleus and part of the redshifted conical outflow.
The projected continuity (in morphology, though not velocity) between the
outflowing gas 1.5$\arcsec$–4$\arcsec$ to the N and the larger scale filaments
(filaments W1 and W2 in the WFM) which are 10s of arcseconds in extent, but
appear to cross the nucleus, supports a scenario in which filament W1 and W2
in the WFM are also in front of the nucleus. With respect to the inner ionized
gas disk: if the M87 jet is intrinsically two-sided, then the (brighter due to
Doppler boosting) NW jet is oriented towards us. For the inner ($<$0.1″)
accretion disk to be perpendicular to jet axis, we expect that the SE half of
the disk is the near side, and the NW 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 SE side is the near side.

Many of the filaments seen in the inner 4″ appear to cross the nucleus. We see
all of the following three cases (a) those which follow the potential of the
black hole, and thus likely truly enter the inner 5 pc; (b) those which 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 highly blueshifted filament participates in a velocity-offset cuasi-
Keplerian rotation as it crosses the nucleus in projection. A potential
explanation for this is gas in a high inclination orbit around the black hole
which is swept up in the biconical outflow.

The twisted isophotes seen in the inner arcsec of the ionized gas disk
potentially resolves 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. 8
provide clear evidence for twisted velocity isophotes in the inner arcsecond
of the disk. The twist is clearest from PA $\sim$42″ at
0$\aas@@fstack{\prime\prime}$3–0$\aas@@fstack{\prime\prime}$5 to PA $\sim$90″
at 0$\aas@@fstack{\prime\prime}$1–0$\aas@@fstack{\prime\prime}$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. 8). Here, the residual over-subtraction is in PA
$\sim$17∘, perpendicular to the jet axis. Since our innermost resolved region
is about $\sim$5 pc, there is of course room for further twists down to the
jet launching region around the black hole. For radii
$\geq$0$\aas@@fstack{\prime\prime}$3, the redshifted side of the disk shows
some indications for further twists in the velocity isophotes, though at 10s
of arcsecond scales we again support a PA close to 40 degrees and inclination
equal to the inner arcsec 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$\alpha$+[N  II] and [O  I] moment 0 maps do not show a clear and
consistent inclined disk shape. In these, while the different intensity levels
are highly non axi-symmetric, 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 N ionized gas in the outflow is seen even in the inner arcsec,
elongating the N-S 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$\aas@@fstack{\prime\prime}$2 to 0$\aas@@fstack{\prime\prime}$6 is slightly
better fit with an inclination of $\sim$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 a i = 42∘ disk). Together, the three results favor a
relatively face-on disk rather than one at i $\sim$ 42∘. This lower
inclination disk is also consistent with the posited 3-dimensional direction
of the jet. Though, 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
which cross the nucleus (sometimes only in projection), and 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 subarcsec-scale disk, and thus black hole mass.

Nevertheless, the deep and integral spectroscopy from MUSE allow us to clarify
and constrain previous competing black hole mass measurements: e.g., G11, W13,
EHT Collaboration, and Liepold et al. (2023). First and foremost, we have
shown that the complexities of the nuclear ionized gas highly complicates the
determination an accurate measurement of the black hole mass from ionized gas
kinematics, which immediately lends more weight to the values measured by G11,
the EHT Collaboration, and Liepold et al. (2023). However, some indications
can be leveraged from the new data, especially via the [O  I] velocities. With
the caveat of all the complexities involved, the model which best fits the [O
I] velocity map and pv diagrams is a G11 or EHT mass black hole together with
a smaller inclination ($\sim$25∘ instead of, e.g., 42∘), i.e. our 6.6i25
model. Comparing the 3.3i42 model and the 6.6i25 model with the observed
velocity data, in most cases we see no strong reason to favor one against the
other as their predictions are very similar within the complexities, and large
linewidths, of the emission lines. However three results support a lower
inclination disk: (a) The [O  I] velocity field residuals (weakly) more
symmetric and less scattered for the 6.6i25 model as compared to the 3.5i42
model; (b) morphologically, the spiral arms in the inner arcsec 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; (c)
the outer ionized gas filaments on 10s of arcsec scales are better fit with a
rotating disk with i = 25∘, rather than i = 42∘. At these large nuclear
distances the rotation is dominated by the galaxy potential, so that there is
no longer a degeneracy between black hole mass and inclination in our
rotational fits. Finally, and independently, we note that the lower
inclination disk provides a more consistent alignment between the axes of the
ionized gas disk, the jet, and the biconical outflow. Thus the most consistent
explanation for the observed inner ionized gas rotating disk is a black hole
mass with mass closer to 6.6 $\times$ 109 surrounded by a low inclination
disk.

To reduce the tension between the G11 and W13 mass measurements, Jeter et al.
(2019) have offered intrinsic sub-Keplerian disk rotation of the ionized gas
as a potential solution. While this remains a viable explanation, especially
for gas in some of the filaments crossing the nucleus, our velocity and
velocity residual results suggest that a lower inclination ionized gas disk is
a simpler explanation for the tension.

We have presented strong evidence for a wide angle (half-opening angle of
$\sim$45∘) outflow along the same (3 dimensional) axis as the jet. This is
supported by the velocity and velocity gradients seen in the blueshifted (N)
and redshifted (S) ionized gas seen in the regions 1$\arcsec$–4$\arcsec$ from
the nucleus. With a radial velocity of $\sim$400 km s-1, this requires to be
powered by the AGN, but since the blue and red shifted cones are only
partially filled with ionized gas, the ionized gas would have to be entrained
from the patchy and filamentary ISM by a wide angle wind launched near the
black hole. Each cone requires to be at least 3″ in extent, and thus the
crossing time is 0.6 Myr, significantly shorter than the lifetime of the M87
larger scale jet and lobes (Owen et al., 2000).

Putting our results together gives us a simple and consistent picture of the
nucleus of M87. A black hole mass closer to $\sim 6.6$ (rather than 3.5)
billion solar masses (for a distance of 17.9 Mpc), together with a warped (in
PA at least) 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 line of sight 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$\aas@@fstack{\prime\prime}$1 to
30″. Schwarzschild modelling of the observed intensity, velocity, and line of
sight velocity distribution of the stars is deferred to Osorno et al., in
preparation. Further, 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 pc-scale to kpc-
scale radio jet. These are being explored in Richtler et al., in preparation.

###### Acknowledgements.

This work was funded by the National Agency for Research and Development
(ANID) / Scholarship Program / Doctorado Nacional / 2017-21171690. We
acknowledge funding from ANID Chile via Nucleo Milenio TITANs (NCN19$-$058),
Fondecyt 1221421, and Basal FB210003. Based on observations collected at the
European Southern Observatory under ESO programme 0103.B-0581.

## References

  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, 7735, 773508
  * Barth et al. (2016) Barth, A. J., Darling, J., Baker, A. J., et al. 2016, ApJ, 823, 51
  * Bentz et al. (2010) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2010, ApJ, 716, 993
  * Bittner et al. (2019) Bittner, A., Falcón-Barroso, J., Nedelchev, B., et al. 2019, A&A, 628, A117
  * Boselli et al. (2019) Boselli, A., Fossati, M., Longobardi, A., et al. 2019, A&A, 623, A52
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798
  * Cappellari & Copin (2003) Cappellari, M. & Copin, Y. 2003, MNRAS, 342, 345
  * Cappellari & Emsellem (2004) Cappellari, M. & Emsellem, E. 2004, PASP, 116, 138
  * Davis et al. (2013) Davis, T. A., Alatalo, K., Bureau, M., et al. 2013, MNRAS, 429, 534
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
  * Emsellem et al. (2004) Emsellem, E., Cappellari, M., Peletier, R. F., et al. 2004, MNRAS, 352, 721
  * Emsellem et al. (2014) Emsellem, E., Krajnovic, D., & Sarzi, M. 2014, MNRAS, 445, L79
  * Event Horizon Telescope Collaboration (2019) Event Horizon Telescope Collaboration. 2019, ApJ, 875, L6
  * Falcón-Barroso et al. (2006) Falcón-Barroso, J., Bacon, R., Bureau, M., et al. 2006, MNRAS, 369, 529
  * Ferrarese & Merritt (2000) Ferrarese, L. & Merritt, D. 2000, ApJ, 539, L9
  * Ford et al. (1994) Ford, H. C., Harms, R. J., Tsvetanov, Z. I., et al. 1994, ApJ, 435, L27
  * Gavazzi et al. (2000) Gavazzi, G., Boselli, A., Vílchez, J. M., Iglesias-Paramo, J., & Bonfanti, C. 2000, A&A, 361, 1
  * Gebhardt et al. (2011) Gebhardt, K., Adams, J., Richstone, D., et al. 2011, 217, 422.05
  * Gebhardt & Thomas (2009) Gebhardt, K. & Thomas, J. 2009, ApJ, 700, 1690
  * Gralla et al. (2020) Gralla, S. E., Lupsasca, A., & Marrone, D. P. 2020, Phys. Rev. D, 102, 124004
  * Gültekin et al. (2009) Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198
  * Harms et al. (1994) Harms, R. J., Ford, H. C., Tsvetanov, Z. I., et al. 1994, ApJ, 435, L35
  * Jeter & Broderick (2021) Jeter, B. & Broderick, A. E. 2021, ApJ, 908, 139
  * Jeter et al. (2019) Jeter, B., Broderick, A. E., & McNamara, B. R. 2019, ApJ, 882, 82
  * Johannsen & Psaltis (2010) Johannsen, T. & Psaltis, D. 2010, ApJ, 718, 446
  * Kormendy & Gebhardt (2001) Kormendy, J. & Gebhardt, K. 2001, in American Institute of Physics Conference Series, Vol. 586, 20th Texas Symposium on relativistic astrophysics, ed. J. C. Wheeler & H. Martel, 363–381
  * Liepold et al. (2023) Liepold, E. R., Ma, C.-P., & Walsh, J. L. 2023, Keck Integral-Field Spectroscopy of M87 Reveals an Intrinsically Triaxial Galaxy and a Revised Black Hole Mass
  * Macchetto et al. (1997) Macchetto, F., Marconi, A., Axon, D. J., et al. 1997, ApJ, 489, 579
  * Marconi & Hunt (2003) Marconi, A. & Hunt, L. K. 2003, ApJ, 589, L21
  * Murphy et al. (2011) Murphy, J. D., Gebhardt, K., & Adams, J. J. 2011, ApJ, 729, 129
  * Owen et al. (2000) Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611
  * Rusli et al. (2011) Rusli, S. P., Thomas, J., Erwin, P., et al. 2011, MNRAS, 410, 1223
  * Saglia et al. (2016) Saglia, R. P., Opitsch, M., Erwin, P., et al. 2016, ApJ, 818, 47
  * Sargent et al. (1978) Sargent, W. L. W., Young, P. J., Boksenberg, A., et al. 1978, ApJ, 221, 731
  * Sarzi et al. (2006) Sarzi, M., Falcón-Barroso, J., Davies, R. L., et al. 2006, MNRAS, 366, 1151
  * Sarzi et al. (2018) Sarzi, M., Spiniello, C., La Barbera, F., Krajnović, D., & van den Bosch, R. 2018, MNRAS, 478, 4084
  * Schnorr Müller et al. (2011) Schnorr Müller, A., Storchi-Bergmann, T., Riffel, R. A., et al. 2011, MNRAS, 413, 149
  * Smith et al. (2021) Smith, M. D., Bureau, M., Davis, T. A., et al. 2021, MNRAS, 503, 5984
  * Vazdekis et al. (2010) Vazdekis, A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2010, MNRAS, 404, 1639
  * Walker et al. (2018) Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128
  * Walsh et al. (2013) Walsh, J. L., Barth, A. J., Ho, L. C., & Sarzi, M. 2013, ApJ, 770, 86
  * Werner et al. (2010) Werner, N., Simionescu, A., Million, E. T., et al. 2010, MNRAS, 407, 2063
  * Young et al. (1978) Young, P. J., Westphal, J. A., Kristian, J., Wilson, C. P., & Landauer, F. P. 1978, ApJ, 221, 721
  * Zhao et al. (2009) Zhao, J.-H., Morris, M. R., Goss, W. M., & An, T. 2009, ApJ, 699, 186

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

Figure 14: A comparison of the pv diagrams obtained from previous HST/STIS
long-slit spectra and our VLT/MUSE-NFM-AO datacube for the [N  II]
$\lambda$6583 (left two panels) and [O  I] $\lambda$6300 (right two panels)
emission lines. The pv diagrams are along a slit through the nucleus of M87 in
PA = 51∘. The horizontal white lines correspond the recessional velocity of
M87 adopted by W13: Vrecc = 1335 km s-1.

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 factor $\sim$11 larger than HST. In adaptive optics mode, MUSE delivers a
higher spatial resolution than STIS. The STIS data used in Walsh et al. (2013)
consists 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$\arcsec\times$8$\arcsec$. 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 dataset 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 14 compares the pv diagrams obtained from the STIS data of Walsh et al.
(2013) and our MUSE data along a nuclear slit in PA 51∘ (the slit PA used by
W13). Note the improvement of the quality of the data, especially in the
weaker [O  I] $\lambda$6300 line.

##  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 the penalized Pixel
Fitting Technique (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 which 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 SNR in each bin of
300 for the WFM datacube and 80 for the NFMdatacube. The continuum SNR for the
Voronoi binning was calculated within the wavelength range 5700 – 5750 $\AA$,
and spaxels with SNR less than $3$ in the continuum were excluded from the
binning. Bins which 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), that covers 3525 – 7500 $\AA$ and includes the
Mg $\lambda$5177 and Na $\lambda$5896 lines. The GIST Pipeline was run on the
wavelength range 5050 – 6000 $\AA$ for the WFM datacube, following Sarzi et
al. (2018) and Emsellem et al. (2014); and the wavelength range 5240 – 5770
$\AA$ for the NFM datacube, in order to avoid the Mg $\lambda$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 LOSVD. We used a 6th (10th) degree additive
polynomial for the NFM (WFM) datacube. Note 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 4 moments to be extracted and used a 4th-degree
multiplicative polynomial. The wavelength ranges were the same as used in the
stellar kinematics run.

The ionized gas kinematics can be also modelled 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 $\AA$, since it
covers the main emission lines in the MUSE range.

While binning improves the SNR 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 use the emission line moment maps which 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 $\sim$1000 km s-1,
relatively close line complexes, e.g., H$\alpha$+[N  II]
$\lambda\lambda$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] $\lambda$4959/[O
III] $\lambda$5007 $=0.35$; [O  I] $\lambda$6364/[O  I] $\lambda$6300
$=0.333$; and [N  II] $\lambda$6548/[N  II] $\lambda$6583 $=0.34$.

##  Appendix C PV diagrams of NFM filaments

Figure 15: As in Fig. 5 but for the remaining four of the eight filaments
marked in the NFM moment maps (black dotted curves in Fig. 4).

Figure 15 shows the pv diagrams of the NFM filaments presented in Sect. 5.2.2,
but not shown in Fig. 4.

##  Appendix D Biconical outflow models

Figure 16: As in the middle and right panels of Fig. 13, but this time showing
both cones superimposed in a single panel. The left panel illustrates a
geometry which allows blueshifts to the N and redshifts to the S, as seen in
the left panel of Fig. 13. 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. in PA $-$72∘.

Figure 16 shows two models of a biconical outflow, to illustrate alternative
scenarios to the biconical outflow discussed in Sect. 5.4 and shown in Fig.
13. Here we show the velocity field expected if both cones are fully filled
with ionized gas. As the observed blueshifted region to the N is wider than
the redshifted region to the S, 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 red shifted 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.