# Observations of AGN-driven feedback: dynamics and ionization of the filaments in M87

Camille Poitras, Marie-Lou Gendron-Marsolais, Valeria Olivares, Yuan Li, Adrien Picquenot, Aurora Simionescu, Matteo Fossati, Alessandro Boselli, Laura Hermosa Muñoz, Sara Cazzoli, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Camille Poitras, Marie-Lou Gendron-Marsolais, Valeria Olivares, Yuan Li, Adrien Picquenot, et al. Observations of AGN-driven feedback: dynamics and ionization of the filaments in M87. Monthly Notices of the Royal Astronomical Society, 2026, 545 (2), 10.1093/mnras/staf1981. hal-05477344

HAL Id: hal-05477344

https://hal.science/hal-05477344v1

Submitted on 20 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Observations of AGN-driven feedback: dynamics and ionization of the filaments in M87

Camille Poitras<sup>1,★</sup> Marie-Lou Gendron-Marsolais<sup>1</sup>,<sup>1</sup> Valeria Olivares<sup>1,2,3</sup> Yuan Li<sup>1,4</sup>,<sup>5</sup> Adrien Picquenot<sup>1,5</sup>, Aurora Simionescu<sup>1,6</sup>, Matteo Fossati<sup>1,7</sup>, Alessandro Boselli,<sup>8</sup> Laura Hermosa Muñoz<sup>1,9</sup>, Sara Cazzoli,<sup>10</sup> Julie Hlavacek-Larrondo<sup>11</sup> and Annabelle Richard-Laferrière<sup>12</sup>

Accepted 2025 November 6. Received 2025 November 6; in original form 2025 September 4

## ABSTRACT

We present a comprehensive kinematic and ionization analysis of the warm ionized filaments (~ 10⁴ K) in M87, the central galaxy of the Virgo cluster, using new integral field spectroscopy from MEGARA (Multi-Espectrógrafo en GTC de Alta Resolución para Astronomía, GTC) and SITELLE (Spectromètre Imageur à Transformée de Fourier pour l'Etude en Long et en Large de raies d'Émission, CFHT). MEGARA targets the southeastern (SE) filaments (~3 kpc from the nucleus) coincident with the only known molecular gas clump, and the far eastern (FE) filament (~15 kpc), spatially isolated within an old radio lobe. SITELLE fully maps the filaments, offering the first complete views of their kinematics and excitation. Combined with archival ALMA (Atacama Large Millimeter Array), MUSE (Multi-Unit Spectroscopic Explorer), and Chandra data, these observations offer a multiphase view of gas dynamics. The filaments display complex motions inconsistent with simple rotation. Velocity structure functions (VSFs) of the warm and cold gas in the central and SE filaments show consistent steep slopes (> 2/3) and flattening on small scales of a few hundred parsecs, possibly suggesting energy injection from Type Ia supernovae, though interpretation is method-limited. The FE filament shows lower VSF amplitude, suggesting less active driving. ALMA CO emission is spatial and kinematically aligned with the ionized gas, the latter showing broader velocity dispersions. Ionization diagnostics indicate AGN (active galactic nucleus)-related processes (e.g. shocks) dominate, with higher energy excitation near the radio lobes and (lower energy) fossil feedback signatures in the FE filament. Finally, the filaments follow the same strong Hα–X-ray surface brightness correlation seen in other clusters, supporting local thermal coupling between phases. However, the FE filament deviates from this trend, possibly due to uplift from past AGN outbursts or limitations in the analysis method.

**Key words:** galaxies: active – galaxies: individual: M87 – galaxies: kinematics and dynamics.

## 1 INTRODUCTION

In cool-core galaxy clusters, where the central cooling time of the intracluster medium (ICM) is shorter than the Hubble time, the hot X-ray emitting plasma is expected to cool and condense. However, observed cold gas masses and star formation rates fall well below classical cooling flow predictions (see e.g. J. R. Peterson & A. C.

Fabian 2006 for a review), although some studies suggest that significant cooling may remain hidden from direct detection (e.g. A. C. Fabian et al. 2022). This discrepancy is now attributed to active galactic nucleus (AGN)-driven feedback, where mechanical and radiative processes heat the ICM, limiting condensation into a multiphase medium (e.g. B. R. McNamara & P. E. J. Nulsen 2007; J. Hlavacek-Larrondo, Y. Li & E. Churazov 2022 for reviews; M. Gaspari, P. Temi & F. Brighenti 2017). Nevertheless, how AGN energy is transported, dissipated, and coupled to different gas phases remains unclear.

<footnote n="1"><sup>1</sup>Département de Physique, de Génie Physique et d'Optique, Université Laval, Québec (QC), G1V 0A6, Canada</footnote>

<footnote n="2"><sup>2</sup>Departamento de Física, Universidad de Santiago de Chile, Santiago 9170124, Chile</footnote>

<footnote n="3"><sup>3</sup>Center for Interdisciplinary Research in Astrophysics and Space Exploration (CIRAS), Universidad de Santiago de Chile, Santiago 9170124, Chile</footnote>

<footnote n="4"><sup>4</sup>Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>

<footnote n="5"><sup>5</sup>Department of Astronomy, Louisiana State University, Baton Rouge, LA 70803, USA</footnote>

<footnote n="6"><sup>6</sup>SRON Space Research Organization Netherlands, Niels Bohrweg 4, NL-2333CA Leiden, The Netherlands</footnote>

<footnote n="7"><sup>7</sup>Dipartimento di Fisica G. Occhialini, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="8"><sup>8</sup>Aix-Marseille Université, CNRS, CNES, LAM, Marseille F-13388, France</footnote>

<footnote n="9"><sup>9</sup>Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Madrid, Spain</footnote>

<footnote n="10"><sup>10</sup>Instituto de Astrofísica de Andalucía, IAA-CSIC, Glorieta de la Astronomía S/N, Granada E-18008, Spain</footnote>

<footnote n="11"><sup>11</sup>Département de Physique, Université de Montréal, Succ. Centre-Ville, Montréal, Québec, H3C 3J7, Canada</footnote>

<footnote n="12"><sup>12</sup>Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote>*E-mail: camille.poitras.2@ulaval.ca</footnote>
## 2 C. Poitras et al.

Warm ionized filaments (~$10^4$ K) are common in brightest cluster galaxies (BCGs) and are believed to trace the dynamic interplay between AGN feedback and the ICM (e.g. M. Gaspari et al. 2018; Y. Li et al. 2020; V. Olivares et al. 2023). Extending tens of kiloparsecs, these structures are typically spatially with soft X-ray-emitting plasma (~$10^7–10^8$ K) and cold molecular gas (~10–100 K), and often surround radio lobes (e.g. P. Salomé et al. 2008; H. R. Russell et al. 2017; V. Olivares et al. 2022a, b, 2025). Despite extensive observations, their origin remains debated. One leading scenario suggests that the filaments condense *in situ* from the ICM via thermal instabilities. In several models, this is triggered by the uplift of low-entropy gas by AGN-inflated bubbles to altitudes where the gas becomes thermally unstable, leading to precipitation that can rain back toward the centre, triggering the supermassive black hole (SMBH; e.g. F. Pizzolato & N. Soker 2005; M. McCourt et al. 2012; M. Gaspari, M. Ruszkowski & S. P. Oh 2013; Y. Li et al. 2015; G. M. Voit et al. 2015; B. R. McNamara et al. 2016; G. M. Voit et al. 2017). Alternatively, external processes such as gas-rich mergers or ram-pressure stripping may deliver multiphase material to the cluster core (W. B. Sparks, H. C. Ford & A. L. Kinney 1993; M. L. Weil, J. Bland-Hawthorn & D. F. Malin 1997; L. Mayer et al. 2006).

M87 (Virgo A, NGC 4486, 3C274), the central galaxy of the nearby Virgo Cluster (~16.5 Mpc; M. Cantiello et al. 2024), provides a well-resolved environment for studying these processes. It hosts a ~$6.5 \times 10^9$ M$_\odot$ SMBH (Event Horizon Telescope Collaboration 2019), powerful relativistic jets with extended radio lobes (D. C. Hines, F. N. Owen & J. A. Eilek 1989; F. N. Owen, J. A. Eilek & N. E. Kassim 2000), and complex X-ray structures, including elongated hot gas filaments within a diffuse halo (e.g. A. J. Young, A. S. Wilson & C. G. Mundell 2002; W. Forman et al. 2007; E. Churazov et al. 2008; N. Werner et al. 2010). M87 is also surrounded by a prominent multiphase filamentary nebula first identified decades ago (H. C. Arp 1967; H. C. Ford & H. Butcher 1979; W. B. Sparks et al. 1993). The warm ionized component extends out to ~18 kpc from the nucleus, exhibiting complex morphologies that closely trace buoyantly rising X-ray cavities and radio structures (A. Boselli et al. 2019, hereafter B + 19), with spatially ultraviolet (UV) filaments of similar extent (W. B. Sparks et al. 2009; P. Tamhane et al. 2025). Notably, M87 is unusually deficient in molecular gas, hosting only a compact CO-bright clump ~3 kpc southeast of the nucleus (A. Simionescu et al. 2018). Moreover, some filaments appear detached from the nebulae (G. Gavazzi et al. 2000, hereafter G + 00, B + 19) and current AGN activity. These may preserve signatures of past feedback events and offer unique insight into the long-term impacts of AGN-driven outbursts.

Despite a rich history of multiwavelength studies (e.g. W. B. Sparks et al. 1993; E. Churazov et al. 2001; A. J. Young et al. 2002; W. B. Sparks et al. 2004; W. Forman et al. 2007; E. Churazov et al. 2008; N. Werner et al. 2013; W. Forman et al. 2017; A. Simionescu et al. 2018; A. Boselli et al. 2019), a comprehensive dynamical characterization of M87's warm ionized filament network remains incomplete. As a result, fundamental questions persist regarding filament formation, motion, and mechanisms by which they trace or respond to both current and past AGN feedback.

In this context, we present new integral field spectroscopy of two key regions in M87's filamentary nebula with MEGARA (Multi-Espectrógrafo en GTC de Alta Resolución para Astronomía). The first targets the inner region (~3 kpc southeast of the nucleus) beyond the inner radio lobes and overleaping the only known molecular gas clump. This region lies inside and near the edge of a shock cocoon propagating through M87's X-ray atmosphere (E. T. Million et al.

2010; A. Simionescu et al. 2018). The second probes the outer eastern filaments (~15 kpc from the nucleus), spatially isolated within an old lobe. Both regions fall outside the previous MUSE wide-field mode (WFM) pointing (1 arcmin × 1 arcmin). Complementary SITELLE (Spectromètre Imageur à Transformée de Fourier pour l'Etude en Long et en Large de raies d'Emission) imaging Fourier transform spectroscopy provides full mapping of the ionized gas across the network. Our aims are to: (1) characterize ionized gas kinematics and turbulence via velocity structure functions (VSFs) and assess coupling with molecular gas; (2) perform emission-line diagnostics; and (3) revisit the Hα–X-ray spatial correlation at large radii, building on V. Olivares et al. (2025).

This paper is structured as follows: Section 2 describes MEGARA and SITELLE observations, data reduction, and archival data sets. Results addressing our three main goals are presented and discussed in Sections 3–5, with conclusions in Section 6. For consistency with previous studies (N. Werner et al. 2013; A. Simionescu et al. 2018), we adopt a systemic velocity of $v = 1307$ km s⁻<fnref n="1" /> for M87 ($z = 0.00436$, $D_L = 16.1$ Mpc; J. P. Blakeslee et al. 2009), where 1 arcmin = 4.68 kpc.

## 2 OBSERVATIONS AND DATA REDUCTION

This study is based on a comprehensive data set, incorporating both new observations and archival data used for the core analysis. A summary of these data is provided in Table 1. Details of the new data and their reduction are described in Sections 2.1 and 2.2, while archival data are documented in the cited references. An overview of M87, its warm ionized filaments, and the main instrument pointings used in this study is shown in Fig. 1.

### 2.1 MEGARA

The observations were conducted in 2024 February (PI: M. Gendron-Marsolais) using the optical integral field unit (IFU) MEGARA in LCB (large compact bundle) mode. MEGARA is installed on the Gran Telescopio Canarias (GTC), a 10.4 m telescope at the Roque de los Muchachos Observatory. The IFU consists of 567 fibres (each with a diameter of 100 µm) arranged in a square microlens array, covering a field of view (FOV) of 12.5 arcsec × 11.3 arcsec. Each microlens has a hexagonal shape inscribed within a projected 0.62 arcsec circular diameter. For these observations, LR-R Volume-Phased Holographic (VPH, i.e. VPH675-LR) was used, covering 6096.5–7303.2 Å with a spectral resolution of $R \sim 6000$, corresponding to a velocity resolution of $\Delta v \sim 50$ km s⁻<fnref n="1" />. The average seeing was 1.0 arcsec (~78 pc).

Two regions in M87 were targeted, sampling different ionized-gas environments. The southeastern filaments (SE filaments, ~3 kpc from the nucleus) form a knotty complex and host the only known detection of molecular gas in M87 (A. Simionescu et al. 2018). The far eastern filament (FE filament, ~15 kpc), the 'plume' first identified by G + 00, is completely isolated from the main nebula. The SE filaments were observed in three observing blocks during the same night, while the FE filament was observed in two. The FE pointing was rotated by 60° northeast to match the filament's orientation.

Each observing block was reduced independently using the MEGARA Data Reduction Pipeline MEGARA DRP<fnref n="1" /> (version 0.15). The pipeline performs bias subtraction, flat-field correction, spectra

<footnote n="1">¹https://github.com/guaix-ucm/megaradrp</footnote>
<table><caption>**Table 1.** Summary of M87 observations used in our analysis.</caption><thead><tr><th>Facility| Instrument (mode/filter)</th><th>Exposure time</th><th>Prog./Obs. ID</th><th>References<sup>a</sup></th></tr></thead><tbody><tr><td>CO(2-1)</td><td></td><td></td><td></td></tr><tr><td>ALMA (Band 6/211-275 GHz)</td><td>5.4 h</td><td>2013.1.00862.S</td><td>A. Simionescu et al. (2018)</td></tr><tr><td>[C II]µm</td><td></td><td></td><td></td></tr><tr><td>Herschel Space Observatory| PACS</td><td>10 880 s</td><td>1342000236278</td><td>N. Werner et al. (2013)</td></tr><tr><td>Optical</td><td></td><td></td><td></td></tr><tr><td>GTC| MEGARA (LCB, VPH675-LR)</td><td></td><td>GTCMULTIPLE2H-23B</td><td>*</td></tr><tr><td>SE filaments</td><td>~ 2 h</td><td></td><td></td></tr><tr><td>FE filament</td><td>~ 5 h</td><td>...</td><td>*</td></tr><tr><td>CFHT| SITELLE</td><td></td><td></td><td></td></tr><tr><td>(SN1)</td><td>~ 1.84 h</td><td>2388234</td><td>*</td></tr><tr><td>(SN3)</td><td>~ 5.67 h</td><td>2274329</td><td>*</td></tr><tr><td>VLT| MUSE (WFM)</td><td>2 × 0.5 h</td><td>60.A-9312</td><td>A. Boselli et al. (2019)</td></tr><tr><td>X-ray (0.5-2.0 keV)</td><td></td><td></td><td></td></tr><tr><td>Chandra</td><td>1745 ks</td><td>See V. Olivares et al. (2025), supporting information</td><td>E.g. A. J. Young et al. (2002); W. Forman et al. (2005, 2007)</td></tr></tbody></table>

Notes. * This paper.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Composite image of M87 combining SDSS optical data (Sloan Digital Sky Survey; Abdurro'uf et al. 2022) with VESTIGE Hα + [N II] emission (Virgo Environmental Survey Tracing Ionised Gas Emission; in red), obtained with MegaCam/CFHT (A. Boselli et al. 2018, 2019). The image reveals the central filamentary network and two isolated plumes extending northeastward at projected distances of ~15 and ~18 kpc from the galaxy centre. Previous MUSE WFM observations are indicated by a pink dashed outline. The MEGARA pointings toward the SE and FE filaments are shown in solid cyan boxes. The blue contour highlights the only known CO(2-1) molecular gas, detected in the SE filaments. Dotted grey contours indicate the 90 cm VLA radio continuum (Very Large Array; F. N. Owen et al. 2000). The field correspond to the full width and 60 per cent of the height of the SITELLE FOV, centred on the galaxy.</figcaption></figure>

tracing and extraction, fibre and pixel transmission correction, as well as wavelength and flux calibration.

During observations, MEGARA simultaneously uses 56 dedicated sky-subtraction fibres arranged in eight groups located 1.75–2.0 arcmin from the FOV centre. While the MEGARA DRP

performs sky subtraction, significant residuals remained after processing. These arise because sky spectra vary with fibre positions in the pseudo-slits, which correspond to their locations at the spectrograph entrance. To mitigate these residuals, we implemented an alternative sky subtraction method. Instead of subtracting a

<footnote>a Reference(s) to the paper in which the mentioned data used in this study – or part of it, in the case of X-ray – were first published.</footnote>
single mean sky spectrum from all fibres, we subtracted the mean of the spectra from the two nearest sky fibre groups for each fibre within the pseudo-slit. This approach significantly improved the sky subtraction quality, leaving only minor residuals generally outside the emission lines of interest. We also found pixel-to-pixel shifts up to 0.71 Å in the central position of the [OI]λ6300 sky line, likely due to inadequate wavelength calibration. Following A. Arroyo-Polonio et al. (2024), we treated the [OI]λ6300 sky line as a reference: for each spectrum, we fit it with a single Gaussian to determine its observed position, then subtracted the offset from the entire spectrum. This approach effectively minimizes wavelength calibration uncertainties, ensuring accurate kinematics.

Finally, each observing block was resampled onto a regularized grid with 0.31 arcsec square spaxels, producing data cubes with dimension of 43 × 40 × 4300, corresponding to a total of 1720 spectra per cube. The final two data cubes were obtained by combining the individual cubes from their respective observing blocks, using a median and ensuring the proper alignment of their central positions.

### 2.1.1 Stellar continuum modelling and subtraction

Accurate stellar continuum subtraction is essential to isolate ionized gas emission, especially in early-type galaxies like M87 where stellar absorption dominates the central regions. Although the SE filaments lie at a moderate distance from the centre, measurable stellar absorption features (e.g. Hα and Fe Iλ6495) remain and must be corrected. In contrast, the FE filament exhibits little to no stellar absorption, so continuum subtraction was applied only to the SE data cube.

To measure and remove the stellar continuum, we constructed a representative median spectrum by averaging all spaxels within the FOV, excluding edge regions prone to instrumental artifacts. This was justified by the uniform continuum shape and consistent depth of the stellar absorption lines across the field, in particular the isolated Fe Iλ6495 line free from nebular emission. This approach assumes negligible stellar age or metallicity gradients in the modelled region. The median spectrum was fitted using PPXF algorithm (Penalized Pixel-Fitting; M. Cappellari & E. Emsellem 2004). During the fit, we masked emission lines, strong sky residuals, and telluric bands to avoid contamination. We used MILES stellar templates (Medium-resolution Isaac Newton Telescope Library of Empirical Spectra; A. Vazdekis et al. 2010), based on Padova + 00 evolutionary tracks (J. Falcón-Barroso et al. 2011; P. Sánchez-Blázquez et al. 2014), spanning 3525–7500 Å at 2.5 Å full width at half maximum resolution. Although MEGARA provides a higher spectral resolution ($R \sim 6000$), this mismatch is not critical in our case, as the intrinsic stellar velocity dispersion in M87 ($\sim 250$ km s⁻¹ between 2–12 kpc, decreasing by only $\sim 10$ per cent at larger radii; E. R. Liepold, C.-P. Ma & J. L. Walsh 2023) broadens the absorption lines well beyond the template resolution. The best-fitting PPXF model was interpolated onto the spectral grid of each spaxel. The resulting fit to the median spectrum is shown in Appendix A. To recover the absolute continuum level, the model was rescaled individually for each spaxel to match its mean observed flux. The final data cube was finally corrected by subtracting the resulting stellar absorption spectrum from each spaxel.

### 2.1.2 Emission-line fitting

For each spectrum, the Hα, [N II]λ6548, 6583, and [S II]λ6716, 6731 lines were fitted simultaneously with Gaussian profiles using the LMFIT routine, based on the Levenberg–Marquardt least-squares algorithm (M. Newville et al. 2016). The fits were constrained to share a common velocity shift and line broadening. The [OI]λ6300,

6364 lines were excluded from the fit, as they are too faint to be measured reliably in individual spaxels. None the less, in integrated spectra combining spaxels with Hα signal-to-noise ratio (SNR) greater than 3, [OI]λ6300 emission is more clearly detected, with surface brightness of $\sim 2.9 \times 10^{-17}$ and $\sim 5.5 \times 10^{-17}$ erg s⁻¹ cm⁻² arcsec⁻² for the SE and FE filaments, respectively. In both integrated spectra, the line's velocity and velocity dispersion are consistent with those of the other emission lines.

To fit the emission lines, we allowed for the possibility of multiple kinematic components. A preliminary visual inspection revealed that parts of the SE filaments contain two kinematically distinct components with similar velocity dispersions. To optimize the fitting process and determine the appropriate number of components, a progressive method was employed. Initially, all spectra were fitted with a single component. To assess whether an additional component was required, we applied the statistical method proposed by S. Cazzoli et al. (2018). This method involves measuring the standard deviation of a line-free continuum region $\epsilon_c$ (in this case, between [N II]λ6583 and [S II]λ6716), and comparing it to the standard deviation of the residuals between the fitted and observed emission lines $\epsilon_{line}$. If $\epsilon_{line} > 3 \times \epsilon_c$ (following a 3σ criterion), the addition of a second component is considered. A third component was also tested but never required. A visual inspection of the spectra was also performed to validate the statistical selection of components. Spectra exhibiting two distinct components are primarily located in the northwestern region of the SE filaments' FOV (see Fig. 2). Additionally, a slight asymmetry in the optical emission-line profiles – located near the eastern edge of the CO(2–1) detection – hints at a second component there, though this is uncertain: Gaussian fitting becomes ambiguous when components strongly overlap, producing multiple possible solutions. As a result, two-component fits were only retained when statistically justified, exclusively in the northwestern region. In the surrounding areas, transitions between one and two components appear slightly abrupt, an effect likely reflecting the intrinsic limitations of Gaussian fitting.

Finally, the number of components, denoted by *N*, was statistically validated using the Akaike information criterion (AIC; H. Akaike 2011). In most cases, the two-component fits in the northwestern region are justified when the AIC difference, defined as $\Delta$AIC = AIC$_{N+1}$ – AIC$_N$, exceeds 10 (J.-J. Wei, X.-F. Wu & F. Melia 2016).

For the remainder of the SE filaments analysis, we define the first component as the one with velocities between approximately –170 to –75 km s⁻¹, which is dynamically connected to the rest of the filament where only a single component is detected. The second component corresponds to velocities within approximately –75 to –30 km s⁻¹. Throughout this paper, we use the label *1 comp.* to refer to the maps from single-component fits, including those within the two-component region. *1st comp.* denotes the first component within the two-component region and single-component spaxels elsewhere. *2nd comp.* refers exclusively to the second component in the two-component region.

### 2.1.3 Map generation

For each fitted emission line, we measured the central wavelength, velocity dispersion, and flux, along with their uncertainties. The uncertainties from the fitting procedure, which generally exceed those from instrumental calibration, were fully accounted for. These parameters were used to produce flux, velocity, and velocity dispersion maps. Fluxes were corrected for Galactic extinction assuming
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) MEGARA emission maps of [N II]λ6583 + Hα for the first and second kinematic components of the SE filaments (continuum-subtracted), with north at the top and east to the left, as well as for the FE filament. The blue-shaded region in the northwest outlines the area where the presence of two components if confirmed. Only spaxels with SNR<sub>Hα</sub> > 3 are shown. CO(2–1) flux contours from ALMA (above the 3σ level) are overlaid on the SE filament maps, ranging from ~1 to ~6 mJy beam⁻¹. All scale bars correspond to 5 arcsec. (b) Representative MEGARA spectra from selected spaxels: from left to right, one from the SE filaments with a single component, one with two components, and one from the FE filament. A broken axis has been applied to better visualize the main spectral lines, which are marked by red dotted lines. For each spectrum, the fitted line profile (red solid line) is shown, along with individual components where applicable (blue dashed lines). The bottom panel displays the residuals of the fits, with the shaded region indicating the ±ε<sub>c</sub> interval.</figcaption></figure>

$A_V = 0.063$ mag (NASA/IPAC Extragalactic Database<fnref n="2" />) and $R_V = 3.1$, following the J. A. Cardelli, G. C. Clayton & J. S. Mathis (1989) extinction law, with a wavelength-dependent correction applied to each line. Fig. 2 shows the [N II] + Hα<fnref n="3" /> maps along with examples of fitted spectra. Radial velocities were corrected for barycentric motion. Since the observed velocity dispersions $\sigma_{\text{obs}}$ were consistently broader than the instrumental width $\sigma_{\text{ins}}$, estimated from the width of the [O I]λ6300 sky line, a correction for instrumental broadening was applied using the following relation

$$ \sigma_{\text{line}} = \sqrt{\sigma_{\text{obs}}^2 - \sigma_{\text{ins}}^2}. \qquad (1) $$

## 2.2 SITELLE

M87 was observed with the SITELLE imaging Fourier transform spectrometer, installed on the 3.6-m Canada-France-Hawaii Telescope (CFHT), as a part of two observing programs (PI: Julie Hlavacek-Larrondo, 2018; PI: Annabelle Richard-Laferrière, 2019). SITELLE provides a 11 arcmin × 11 arcmin FOV, covering the

entire filamentary network and its surroundings, with spatial sampling of 0.32 arcsec pixel⁻¹ (L. Drissen et al. 2019). The 2018 observations used the SN3 filter (R ~ 950, 651–685 nm), targeting Hα, [N II]λ6548, 6583, and [S II]λ6716, 6731 emission lines. The typical seeing during these observations was ~ 1.2 arcsec. In 2019, additional observations were obtained using the SN1 filter (R ~ 680, 365–385 nm), which covers the [O II]λλ3727–3729 doublet (unresolved).

Data reduction used SITELLE's dedicated pipelines, ORBS and ORCS<fnref n="4" /> (version 3.0; T. Martin, L. Drissen & S. Prunet 2021). Wavelength calibration was initially carried out with a high-resolution laser source, but instrument flexures during target acquisition can introduce velocity shifts of ~15–20 km s⁻¹ across the FOV, with variations in M87's central region remaining below 5 km s⁻¹. These were corrected by measuring OH sky line central positions within cubes, following T. B. Martin, L. Drissen & A.-L. Melchior (2018). Flux calibration was performed for each filter using stars in the field matched to Gaia Data Release 3 (Gaia Collaboration 2023) sources with calibrated optical spectra. Synthetic photometry was derived by integrating Gaia spectra over the filter transmission curve, then compared to SITELLE fluxes measured from 2D Gaussian fits

<footnote n="2">²https://ned.ipac.caltech.edu/</footnote>

<footnote n="3">³[N II]λ6583 + Hα; this definition will be used throughout the paper.</footnote>

<footnote n="4">⁴https://github.com/thomasorb/orbs</footnote>
of stellar profiles in deep frames (i.e. the sum of all interferogram frames). For the SN3 filter, a linear regression between SITELLE and Gaia-based fluxes yielded a slope of $0.942 \pm 0.092$, indicating that SITELLE fluxes slightly underestimate the true flux. A correction factor of $1/0.942$ was therefore applied. For the SN1 filter, a larger discrepancy was found, with a slope of $0.804 \pm 0.040$ and the corresponding factor was applied accordingly. A detailed description of this type of calibration procedure for SITELLE will be presented in Vicens-Mouret (in preparation).

Prior to the emission-line fitting, sky subtraction was performed using a circular region with a radius of $0.5$ arcmin, centred at $(\alpha: 12h31m08.0636s, \delta: +12h24m59.761s)$. This area contains neither point sources nor detectable source emission and provides a representative sky background comparable to other object-free regions of the FOV.

Stellar continuum subtraction for SITELLE required a spatially adaptive approach due to the noticeable spectral variations across the field. Visual inspection of the SN1 and SN3 data cubes revealed significant changes in the continuum level and absorption line depth with galactocentric radius, while stellar absorption velocities remained uniform in both filters, consistent with M87's low rotation ($\sim 25$ km s⁻¹ beyond 5 kpc, decreasing inward; E. R. Liepold et al. 2023). To account for the spatial variations, the field was divided into concentric annuli centred on the galactic nucleus. Radial boundaries were selected to reflect observed spectral gradients. Within each annulus, spaxels representative of the stellar population were identified based on their low $H\alpha$ emission ($< 5 \times 10^{-19}$ erg s⁻¹ cm⁻²) using a flux map derived from the original data cube. Foreground stars and spaxels including the jet were manually excluded. Each annulus contained at least 600 qualifying spaxels. To construct a representative stellar spectrum for each annulus, the selected spectra were first normalized by their mean flux to prevent bias from bright spaxels, then median-combined. The resulting spectra were fitted with PPXF using MILES templates, fitting the SN1 and SN3 filter simultaneously. The best-fitting models were linearly interpolated to the observed spectral grid of each filter and rescaled to match each pixel's mean flux for the continuum level. Subtracting the appropriately scaled model for each annulus produced the final continuum-corrected cube. An example integrated spectrum and its fitted model are shown in Appendix A.

From the calibrated data cubes, with sky and stellar continuum subtracted, ORCS was used to generate maps of the continuum, emission-line amplitude and flux, radial velocity, and velocity dispersion, along with their associated uncertainties. To enhance the SNR, maps were spatially binned ($2 \times 2$). A barycentric correction was applied to the velocity map. Fluxes were corrected for Galactic extinction following the method described in Section 2.1.3. The $[N\text{II}] + H\alpha$ emission map is shown in Fig. 3. Due to the applied $\text{SNR}_{H\alpha} > 3$ threshold, the faint outermost plume filament first identified by B + 19 is not visible in the map. The same figure also presents the $[\text{OII}]\lambda\lambda3727-3729$ emission map, where faint [O II] emission is detected in several filaments, including the FE filament, consistent with earlier reports by H. C. Arp (1967) and H. C. Ford & H. Butcher (1979). However, spaxel-by-spaxel SNR remains below 3 across most regions. Combined with the unresolved doublet in SN1 data, this prevents the use of [O II] for reliable electron density diagnostics. $H\beta$ and $[O III]\lambda5007$ are outside the SITELLE filter coverage, further limiting ionization diagnostics. Based on the mean SNR of spaxels with $\text{SNR}_{H\alpha} > 3$ and the relation from L. Rousseau-Nepton et al. (2019), the minimum measurable velocity dispersion is $\sim 180$ km s⁻¹ for SN3 and $\sim 375$ km s⁻¹ for SN1. Given these limits, velocity-dispersion results are not used in our analysis, though the map is provided for reference in Appendix B.

## 2.3 Archival data, WCS alignment, and spatial resampling

To compare the kinematic and morphological properties of different gas phases, we incorporate archival data sets with our new observations. We directly compare MEGARA maps of the SE filaments with ALMA CO(2–1) flux, velocity and velocity dispersion maps derived via masked-moment analysis (see A. Simionescu et al. 2018 for details). We also compare indirectly SITELLE velocity field with the MUSE velocity map from B + 19, based on the original data of E. Emsellem, D. Krajnovic & M. Sarzi (2014), corrected to the systemic velocity. The MUSE data have a spaxel size of $0.2$ arcsec, $R \sim 1500-3500$, and a seeing of $0.9$ arcsec. In addition, we directly compare the SITELLE $H\alpha$ map with deep Chandra X-ray data (see Table 1).

Spatial consistency was ensured by aligning the World Coordinate Systems (WCS) for all data sets. A visual inspection revealed slight offsets in the MEGARA and MUSE data relative to the others. These were corrected using instrument-specific alignment methods tailored to each data set and the objectives of the analysis. For MEGARA, which lacks bright point sources in its FOV, alignment was achieved by cross-correlating $H\alpha$ filamentary structures with SITELLE (SNR > 3, unbinned), using the correlate2d function after degrading SITELLE to MEGARA's spaxel grid. The resulting pixel offsets were converted into WCS corrections and applied to MEGARA cubes: $(\alpha: -0.000677, \delta: +0.000135)^\circ$ for the SE pointing and $(\alpha: -0.000364, \delta: +0.0000854)^\circ$ for the FE pointing. For MUSE, alignment was achieved by matching the galactic centre position in its velocity maps to SITELLE, with a final correction of $(\alpha: -0.0001785, \delta: +0.0003881)^\circ$.

To enable direct spaxel-by-spaxel comparisons with ALMA, MEGARA maps were convolved with a Gaussian kernel to match the ALMA beam. The ALMA data reach an rms sensitivity of $0.2$ mJy beam⁻¹ per $50$ km s⁻¹ channel, with a synthesized beam of $1.4$ arcsec $\times 0.7$ arcsec (position angle of $75^\circ$; A. Simionescu et al. 2018). The ALMA data were resampled onto MEGARA's spaxel grid using astropy.reproject function. This enabled the construction of three key comparative maps: (1) the peak-normalized $([N\text{II}] + H\alpha)/\text{CO}$ flux ratio, (2) the line-of-sight (LOS) velocity difference, and (3) the velocity dispersion ratio (see Section 3.2.2).

# 3 FILAMENT KINEMATICS

## 3.1 Velocity field across the filament network

M87's ionized gas filaments form a complex and extended network whose kinematics have been extensively studied for decades. Early long-slit spectroscopy revealed irregular velocities and significant variations across the filaments (T. M. Heckman et al. 1989; W. B. Sparks et al. 1993; N. Werner et al. 2013). In particular, W. B. Sparks et al. (1993) proposed a 3D geometry: the northwestern filaments, near the jet, lie in the foreground and are likely outflowing, while the southeastern filaments, near the counter-jet, are in the background and infalling – a scenario supported by dust morphology and radio lobe orientation.

With the advent of integral field spectroscopy, recent studies have refined this picture. MUSE WFM observations (B + 19) revealed highly disordered velocity fields across the central $\sim 1$ arcmin $\times 1$ arcmin, consistent with turbulence likely driven by AGN feedback from the jet (Y. Li et al. 2020). At smaller scales, J. Osorno et al. (2023, hereafter O + 23) used MUSE narrow-field mode (NFM) data to probe the nuclear region in detail, unveiling complex kinematics. Nevertheless, previous velocity studies lacked
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) SITELLE continuum-subtracted [N II] + Hα emission map showing only spaxels with $SNR_{Hα} > 3$ and velocity error ≤ 30 km s⁻¹. (b) Corresponding [O II]λλ3727–3729 emission map. The main panel displays all spaxels, while the inset shows a zoomed-in view of the central region with only spaxels satisfying $SNR_{[OII]} > 3$. In both panels, the MEGARA FOVs targeting the SE and FE filaments are outlined in blue boxes. The displayed maps represent approximately 10 per cent of SITELLE’s full FOV. All scale bars correspond to 1 arcmin. Maps are shown with 2 × 2 binning and a logarithmic colour scale to enhance contrast between bright and faint emission.</figcaption></figure>

the spatial coverage to trace motion across the full filament system: MUSE WFM covered only the central filaments, while earlier long-slit data extended slightly farther. With SITELLE and MEGARA, we now access both the full spatial extent of the filamentary network and sufficient spatial resolution, including the southeastern and isolated filaments beyond 40 arcsec (~3.1 kpc) from the nucleus. Velocity map derived from SITELLE and velocity/dispersion maps from MEGARA are shown in Figs 4 and 5, respectively.

In the nuclear region, a small number of SITELLE spaxels confirm the complex velocity structure reported by O + 23. Their analysis showed that the so-called ionized gas ‘disc’ is not a simple rotating structure, but a mix of rotating material, filamentary emission with distinct bulk velocities, and features crossing the nucleus in projection – producing irregular velocity gradients exceeding 300 km s⁻¹ on subarcsecond scales. Just beyond the nucleus, extending out to the edge of the MUSE WFM pointing, we find excellent agreement with the LOS velocities presented by B + 19 and earlier long-slit studies, with disordered motions typically reaching up to ±300 km s⁻¹.

In the south, the MEGARA data reveal the detailed velocity field of the SE filaments, including a region with resolved two-component profiles located at the edge of the radio counter-jet. The primary component (1st comp.) has a mean LOS velocity of ~−120 km s⁻¹ and a velocity dispersion of ~55 km s⁻¹. It is kinematically aligned with cold molecular gas traced by CO(2−1) and [C II]λ157µm emission, both previously reported as spatial with the ionized gas (N. Werner et al. 2013; A. Simionescu et al. 2018). A more detailed comparison of these phases is presented in Section 3.2.2. In the FE filament, we measure a mean velocity of ~−60 km s⁻¹, consistent with G + 00, and a low velocity dispersion of ~30 km s⁻¹. The kinematic properties of the SE and FE filaments are further discussed in Section 3.2.1 in the context of velocity structure and turbulence.

Finally, SITELLE allows us to probe the most extended southeastern filaments, which had not been kinematically studied before. These structures appear to skirt buoyant radio bubbles of varying ages

rising through M87’s hot atmosphere (W. Forman et al. 2007). These outer filaments exhibit mostly negative velocities between −250 to 50 km s⁻¹. To assess whether M87’s ionized gas filaments show large-scale rotation, we extend the analysis of O + 23 beyond the MUSE WFM pointing. Their study focused on the central filaments (~1 arcmin × 1 arcmin), comparing observed velocities to analytical thin-disc rotation models. They concluded that while the innermost filaments deviate strongly from rotation, outer filaments (beyond ~10 arcsec) can reasonably be described by rotation in a disc inclined by i ≈ 25°, particularly influenced by the stellar potential. Here, we investigate whether this pattern persists into the more distant southeastern filaments beyond 40 arcsec with SITELLE, which were not included in their analysis. Using the same analytical framework as O + 23, the rotational velocity $v_{\text{rot}}$ at radius r is given by

$$v_{\text{rot}}(r) = \sqrt{\frac{G[M_{\bullet} + M_{\text{bulge}}(r)]}{r}}, \quad (2)$$

where G is the gravitational constant, $M_{\bullet} = 6.0 \times 10^9 \, \text{M}_\odot$ (O + 23) is the mass of the SMBH, and

$$M_{\text{bulge}}(r) = \int_{0}^{r} 4\pi r'^2 \Upsilon(r')v(r')dr', \quad (3)$$

is the enclosed stellar mass within a radius r, with $v(r')$ is the V-band luminosity density profile from K. Gebhardt & J. Thomas (2009), and $\Upsilon(r')$ the stellar M/L ratio fitted by E. R. Liepold et al. (2023). As $v(r')$ is non-analytic, the integral is evaluated numerically. To project the rotational velocity into the LOS, we calculate

$$v_{\text{LOS}} = v_{\text{rot}}(r) \sin i \cos \theta, \quad (4)$$

where i = 25° is the disc inclination, and θ is the azimuthal angle in the disc plane between the major axis (assumed at PA = 45°), adopting the geometric parameters derived by O + 23. We compare this model to the SITELLE velocities via a residual map ($v_{\text{obs}} - v_{\text{model}}$), to highlight spatial regions where the observed ionized gas velocities are
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) SITELLE velocity map showing spaxels with $SNR_{H\alpha} > 3$ and velocity uncertainties $\le 30$ km s$^{-1}$ ($2 \times 2$ binning). (b) Rotating thin-disc model computed for the combined gravitational potential of the central SMBH and the stellar bulge, based on the best-fitting model from O + 23. (c) Residual map showing the difference between the observed velocities and the rotation model. For each panel, faint solid white/black lines centred on the nucleus mark the disc's kinematic P.A. (45°) and its perpendicular axis. MEGARA FOVs targeting the SE and FE filaments are outlined in blue boxes, and the existing MUSE WFM observations is outlined with a white/black dashed box. The colour scale is consistent across all panels and each scale bar corresponds to 1 arcmin.</figcaption></figure>

consistent with a disc-like rotation and where significant deviations occurs. The resulting residual map is shown in Fig. 4. Within the MUSE FOV, residuals match O + 23. Beyond this area, few portions of the filaments follow the general trend expected from disc rotation, with residuals on the order of $\sim \pm 50$ km s$^{-1}$. However, the velocity field is more complex. Notably, the SE and FE filaments, as well as the horseshoe-like filament located at $\sim 90$ arcsec southeast from the nucleus, exhibit substantial deviations ($> 100$ km s$^{-1}$) from model

predictions, suggesting additional local processes (e.g. shocks or dynamical instabilities). In particular, the SE filaments lie just beyond the inner radio lobes, while the FE filament is embedded in an older radio lobe. Additionally, the influence of the biconical outflow identified by O + 23 may also contribute to the disturbed motions in the inner filaments. Therefore, a simple rotating disc model is insufficient to explain the full range of filament motions observed in M87.

## 3.2 Velocity structure functions

VSFs have emerged as a valuable tool for characterizing turbulence in filamentary systems and are increasingly applied across entire FOV to investigate dominant turbulent driving mechanisms. In M87, Y. Li et al. (2020) pioneered the use of VSFs in such optical filaments, using MUSE WFM data of the central region. They identified a driving scale of 1–2 kpc, consistent with the AGN jet and X-ray cavity, linking turbulence to jet-driven activity. Here, we extend this approach by analysing the nature of filament motions using first-order VSFs ($n = 1$). As a two-point correlation function, the VSF is defined as

$$ \langle |v(l)|^n \rangle = \langle \| [v(x + l) - v(x)] \|^n \rangle, \quad (5) $$

which quantifies the average absolute velocity difference between two points, $x$ and $x + l$, separated by a LOS vector $l$. VSFs are typically expressed as a function of the scalar separation $l = |l|$, providing a statistical description of velocity variations across spatial scales. Mean velocity differences $\langle \Delta v \rangle$ are computed by binning separation distances $l$ in logarithmic space: at small scales, they are manually adjusted to account for data discretization, while the larger scales – where the sampling is more continuous – fixed width bins are adopted.

To ensure the reliability of the VSFs, spaxels with velocity uncertainties greater than 20 km s$^{-1}$ were excluded, a threshold above the mean uncertainty that does not significantly affect the results. These higher uncertainty spaxels lie mostly at filament edges where the SNR is lower. Additionally, due to the finite size of the FOV and structures themselves, the number of point pairs decreases with increasing separation. As a result, bins containing fewer than 20 per cent of the maximum number of pairs are marked in grey in the VSF plots to indicate reduced statistical significance. Each VSF plot also includes a vertical marker at half the size of the spatial window used, highlighting scales where interpretations become uncertain, and the seeing scale to illustrate potential atmospheric effects. Further discussion of these limitations and the uncertainties they introduce in our interpretation is provided in Section 3.2.6. For the uncertainties of the mean velocity differences, they were estimated using a combined approach that accounts for both error propagation from individual velocity measurements and statistical variations via bootstrapping.

To guide VSF interpretation, we overlay reference power-law slopes on the plots. Specifically, a slope of $1/3$ is shown as a benchmark for classical Kolmogorov turbulence as expected for an incompressible fluid (A. Kolmogorov 1941). We also include a steeper slope of $2/3$, often matching our measurements and commonly observed in other systems (e.g. Y. Li et al. 2020; S. Ganguly et al. 2023). It is also worth noting that for compressible (supersonic) turbulence, the expected slope is approximately $1/2$ (S. A. Boldyrev 1998). While many VSFs follow a clear power-law scaling at smaller spatial separations, several also exhibit a noticeable flattening at intermediate scales. This flattening indicates the dominant driving scale of turbulence, where turbulent kinetic energy injected at large scales begins transferring to smaller scales
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Velocity dispersion maps and (b) velocity maps of the SE (top) and FE (bottom) filaments, obtained with MEGARA (left) and SITELLE (right), showing only spaxels with $SNR_{H\alpha} > 3$. North is up and east is to the left. For the SE filaments in MEGARA, the main maps display the first component, while separate maps display the second component and, for the velocity map, the single-component fit within the two-component region (outlined). SITELLE maps include only a single velocity component. Colour scales are consistent across maps of the same type. All scale bars correspond to 5 arcsec. (c) VSFs for both pointings, derived from SITELLE and MEGARA data. Shaded areas indicates bins with an insufficient number of pairs ($< 20$ per cent of the maximum pair count). Two reference slopes are shown: 1/3 (Kolmogorov turbulence) and 2/3 (steeper slope). Error bars are displayed, although in some cases they are smaller than the symbol size. The vertical dotted and solid lines indicate the SITELLE seeing limit ($\sim 1.2$ arcsec) and the half the height of the MEGARA FOV, respectively.</figcaption></figure>
before being dissipated by viscosity. We return to this interpretation
in the context of specific cases discussed in the results.

In the following sections, we investigate the VSF derived from
different instruments, substructures, and gas phases (Sections 3.2.1–
3.2.3). A global interpretation is presented in Section 3.2.4, followed
by an application to the turbulent heating rate in Section 3.2.5.
The limitations inherent to the VSF analysis are discussed in
Section 3.2.6.

3.2.1 SE and FE filaments

VSFs were first computed for the SE and FE filaments using MEGARA velocity maps. An initial set of VSFs was derived under the assumption of a single-component velocity fit for both pointings (1 comp.). For the SE filaments, a second VSF was derived from the primary component (1st comp.) of the two-component fit to assess the effect of overlapping velocity structures. A separate VSF for the second component was not computed due to the limited spatial coverage and its large velocity offset relative to the rest of the filament. MEGARA VSFs were then compared to those from SITELLE, restricted to the same FOVs. The results are shown in Fig. 5 (panel c).

All VSFs display a clear power-law behaviour at small scales,
with slopes generally steeper than the Kolmogorov 1/3 value,
approaching ~2/3. In the single-component analysis, both MEGARA
and SITELLE show a flattening at scales of ~0.2-0.3 kpc. The
corresponding velocity amplitudes at the flattening scale are ~30 km
s⁻¹ for the SE filaments and ~15-20 km s⁻¹ for the FE filament.
At larger scales, the VSF continues to rise, although the declining
number of spaxel pairs at large separations limits our ability to draw
firm conclusions in this regime.

Interestingly, for the SE filaments, the VSF based on the first com-
ponent does not exhibit a clear flattening; instead, the power increases
steadily toward larger scales. This suggests that the flattening seen
in the one component VSF may not correspond to a physical scale
break, but rather the impact of LOS complexity. Further discussion
of these effects is provided in Section 3.2.6.

Overall trends are similar between MEGARA and SITELLE,
though SITELLE yields slightly higher VSF amplitudes, which may
reflect the increased noise of the SITELLE velocity maps in this
study (see Section 3.2.6).

3.2.2 Comparison between SE filaments and CO(2−1)

To assess whether a spatial and kinematic correspondence exists
between the SE ionized filaments observed with MEGARA and the
CO(2−1) emission detected by ALMA, we directly compare the
two data sets. For this purpose, we use the smoothed and resampled
flux, velocity, and velocity dispersion maps for both gas phases,
along with the corresponding difference and ratio maps as described
in Section 2.3. Part of the CO(2−1) emission overlaps the region
where MEGARA detected two velocity components. The ALMA
observations indicate an LOS velocity of ~130 km s⁻¹, which falls
within the range of MEGARA’s first component. We therefore use
the first component in this region and the single-component map
elsewhere. All results discussed in this section are presented in Fig. 6.

Spatial correspondence was first assessed via the normalized flux
ratio map of ([N II] + Hα)/CO(2−1). The emission peaks of the
first component of the SE filament and the CO(2−1) aligned closely,
yielding ratios near or slightly below unity. Moving eastward, the
ratio increases substantially, with a median value around 4 and

even higher values at the edges. Outside the CO detection, strong
ionized emission persists without molecular counterpart. This may
reflect the destruction or excitation of molecular gas by AGN-driven
processes, as suggested by A. Simionescu et al. (2018). In contrast,
[C II] emission coincides spatially with the ionized filaments and
also extend beyond the CO detection, as reported by N. Werner et al.
(2013).

Kinematic correspondence was evaluated through LOS velocity-
difference and velocity dispersion ratio maps (Hα<fnref n="5" />/CO) to highlight
spatial regions of agreement or discrepancy. We complemented this
analysis with spaxel-by-spaxel scatter plots comparing velocity and
velocity dispersion between the two tracers, following the approach
of G. R. Tremblay et al. (2018).

The velocity-difference map reveals a strong correspondence
across most of the overlapping region, with offsets typically below
55 km s⁻¹ and a mean of (29 ± 1) km s⁻¹. This supports the idea that
ionized and molecular gas are largely cospatial and kinematically
aligned. A similar trend is observed in [C II] emission, which varies
from approximately −100 km s⁻¹ in the southeast to −120 km s⁻¹
northwest along the SE filaments (N. Werner et al. 2013), roughly
consistent with MEGARA Hα velocities spanning approximately
from −80 to −150 km s⁻¹. For the CO, larger velocity offsets
appear at the eastern and western edges of the CO detection. In
the east, the discrepancy may reflect the limits of a single-component
Gaussian fit, consistent with previous indications of a second velocity
component in this region (see Section 2.1.2). In the west, using two
components reduces the mean offset from approximately 70 to 50 km
s⁻¹, suggesting a more accurate kinematic decomposition. A spaxel-
by-spaxel LOS velocity comparison (see Fig. 6, panel a) confirms
the alignment: most points fall near the 1:1 relation, with some
dispersion toward 2:1. Notably, the use of two velocity components
in MEGARA improves the agreement. Given this consistency, we
further compare turbulent properties using VSFs constructed from
both velocity maps (i.e. MEGARA, and ALMA). On small scales,
CO(2−1) and Hα VSFs follow similar power-law trends with slopes
close to 2/3. When the Hα VSF is restricted to the MEGARA spaxels
within the CO(2−1) detection footprint, the agreement between the
two improves. Within this region, both VSFs show similar power-
law slopes and amplitudes at small scales. Beyond ~0.1 kpc, the CO
VSF steepens, yielding higher amplitude than Hα. In contrast, the Hα
VSF maintains a similar slope to that of the full first-component VSF,
remaining roughly parallel. Around 0.2–0.3 kpc – near the half the
width of the CO-emitting region (~300 pc east–west) – the CO VSF
flattens, while the Hα VSF steepens, and their amplitudes converge
again at larger separations where pair statistics become insufficient.

Despite the overall velocity alignment, the velocity dispersions
differ: Hα profiles are systematically broader than CO(2−1), with
mean values of (57 ± 1) km s⁻¹ versus (26 ± 2) km s⁻¹, respectively.
The average ratio is approximately 2:1, rising to 5:2 toward the
edges – though these extreme values may be affected by low SNR
and should be treated with caution. Spaxel-by-spaxel comparison
confirms this trend (Fig. 6, panel a). While velocity centroids
align closely, dispersion ratios range from 1:1 to 4:1, indicating
consistently broader ionized gas line compared to molecular gas. This
behaviour is consistent with previous studies of multiphase filaments
in BCGs (e.g. G. R. Tremblay et al. 2018; V. Olivares et al. 2019),
which attribute broader Hα lines to the higher volume filling factor
of the warm ionized medium increasing the likelihood of multiple

<footnote n="5">⁵Hα is used here as a proxy for the ionized gas. All emission lines in a given
spectrum share the same kinematic profile in our fits.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Spaxel-by-spaxel comparison of ALMA CO(2–1) and MEGARA $H\alpha$ (1st comp.) LOS velocities (blue, bottom left) and velocity dispersion (pink, upper right), smoothed using a Gaussian kernel. For the velocity comparison, the dotted contours represent the distribution of spaxels using the 1 comp. fit, illustrating how the use of a two-component decomposition in the appropriate region tightens the overall distribution. Shaded bands indicate regions where the $H\alpha$ values are 1–2 times (light grey) and 1–4 times (dark grey) the CO values. The 1:1 relation is highlighted by the solid grey line. (b) VSFs for the full SE filament using MEGARA (blue squares), the MEGARA spaxels within the CO(2–1) detection region (black circles) and the CO(2–1) detection itself (magenta unfilled circles). Two reference slopes are shown: 1/3 (Kolmogorov turbulence) and 2/3 (steeper slope). The vertical dotted, dashed, and solid lines indicate the ALMA beam size (1.4 arcsec, lager beam axis), the MEGARA seeing limit (~1 arcsec) and the half-width of the CO-detected region (extent from east to west), respectively. (c) Velocity difference map ($H\alpha$ – CO), with the area right of the white line corresponding to the two-component region. (d) Velocity dispersion ratio map ($H\alpha$/CO); edges should be disregarded. (e) Normalized ([NII] + $H\alpha$)/CO flux ratio map. All scale bars correspond to 5 arcsec.</figcaption></figure>

LOS velocity components. By comparison, the velocity dispersion measured for [C II] in the southern filaments (~55 km s⁻¹ after instrumental resolution correction; N. Werner et al. 2013) aligns well with the $H\alpha$ values. While spectral resolution differences between instruments may partly contribute to this discrepancy, comparable velocity dispersions for the warm filaments are reported with similar resolution data (B. Vigneron et al. 2024). Finally, as suggested by A. Simionescu et al. (2018), the low velocity dispersion of CO may reflect its embedding in a predominantly laminar flow.

### 3.2.3 Across the central filaments

While VSFs are often applied globally, recent studies highlight their ability to reveal localized differences. For example, Y. Li et al. (2020) and S. Ganguly et al. (2023) have demonstrated that VSFs can also capture distinctions between inner and outer filament

regions, revealing how the nature of turbulence may vary with distance from the galaxy centre. Another factor influencing VSF measurement is the viewing geometry. H. Hu et al. (2022) explored how projection effects – both in simulations and observations – can lead to directional dependencies in VSFs. By analysing VSFs across different viewing angles and azimuthal regions in the plane of the sky, they showed that the VSF slope can appear significantly steeper when the LOS is aligned with the outflow direction, and flatter when viewed perpendicularly. Temporal evolution is another factor: M. Fournier et al. (2025) found that the VSF slope evolves over time, steepening during periods of AGN activity, such as jet-driven cavity inflation.

Given these known sources of variability in VSFs, and the localized signatures observed in the SE/FE filaments (Section 3.2.1), we are motivated to investigate whether small-scale variations are present across other regions of the filamentary network. For this purpose, we use the MUSE velocity map (see Fig. 7, panel a), which
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) MUSE velocity map of the central filaments. Each marker corresponds to the centre of a VSF computed within a circular aperture ($R = 0.39$ kpc), shown as the white aperture, and distributed across subregions F1 to F5 and I1 to I5 (as defined in the reference map, lower right). Crosses indicate locations where no VSF flattening is observed. (b) VLA 6 cm radio emission with MUSE filaments contours overlaid in white. The schematic at bottom right shows the definition of the angle relative to the radio jet axis. Panels (c) and (d) present the VSF parameters – top to bottom, slope (black), amplitude (red), and flattening scale (blue) – as functions of projected radial distance from the nucleus and projected angle with respect to the jet axis, respectively. Semitransparent horizontal lines show the median values of each parameter: slope ~0.92 dex, amplitude ~32 km s⁻¹, and flattening scale ~0.18 kpc (excluding data points within the central grey shaded and magenta-highlighted regions). In panel (c), a break in the horizontal axis separates measurements for the SE and FE filaments (star markers). The vertical dashed green line marks the projected length of the optical jet (~1.5 kpc), while the central grey shaded region represents the 0.78 kpc diameter of the VSF aperture. In panel (d), vertical dashed and dotted lines correspond to angles of 90° and 180°, respectively. Highlighted points in panels (a) through (d) are shown in magenta. Filled symbols indicate regions projected near the jet axis, while hollow symbols identify filament segments with velocities significantly offset from their local environments.</figcaption></figure>

captures the central filaments and the area through which most of the visible jet/counter-jet extends. This data set provides both the spatial coverage and a sufficient number of spaxels needed to compute localized VSFs across multiple subregions, enabling us to probe how turbulence varies locally with galactocentric distance and its relation to physical features such as the jet/counter-jet.

Ideally, studying localized turbulence would require isolating individual filaments, but this is challenging due to projection effects. To reduce contamination from overlapping structures, we identified six filamentary subregions based on their projected morphology and kinematic continuity in the MUSE velocity map. This segmentation aims to separate structures that may be distinct in three dimensions
yet appear close in projection and exhibit different velocities (see Appendix C for a comparison of VSFs with and without this separation). The segmentation is inherently subjective, relying on visual identification of morphological and kinematic continuity in projection. Ambiguity is unavoidable – particularly at junctions where filaments overlap or intersect – making it difficult to define strict boundaries. These uncertainties should be kept in mind when interpreting the resulting VSFs. In the following, we describe the subregions in order of increasing projected distance from the galaxy centre, following the nomenclature of W. B. Sparks et al. (2004), without implying any specific infalling or outflowing motion:

F1. A prominent filament extending from the galaxy centre northwestward before curving east. The structure then forms a knot (W. B. Sparks et al. 2004), from which it diverges into a double-stranded morphology. The portion of the knot sharing similar velocities with the previous part of the filaments has been included in this subregion, as well as the southern strand that continues at an angle of ~230° – aligned with the direction of the previous straight segment. This strand gradually fades beyond the MUSE FOV. Overall, this subregion constitutes the longest coherent structures in projection.

F2. A second filament emerging from the galaxy centre and bending northeastward, forming a loop-like structure that appears to rejoin F1 just north of the nucleus in projection.

F3. A northern feature extending vertically from the centre, composed of few narrow strands, which is more clearly seen in [N II] + Hα images with the Hubble Space Telescope (HST – N. Werner et al. 2010; P. Tamhane et al. 2025) and the ESO New Technology Telescope (F. Macchetto et al. 1996; W. B. Sparks et al. 2004).

F4. A small filament that connects the eastern end of F2 with F1, close to the knot. This structure shows a distinct kinematic signature, with a velocity of ~300 km s⁻¹, offset from neighboring filaments (~100 km s⁻¹).

F5 and F6. Located at the southeastern edge of the MUSE FOV. F6 is associated with the SE filaments. F5 appears kinematically distinct from F1, with a velocity of ~−30 km s⁻¹ compared to ~80 km s⁻¹ in adjacent regions. F5 traces the easternmost strand of the double-stranded morphology emerging after the knot described in F1, and gradually fades outside the MUSE FOV.

We also identified five small, spatially isolated filaments in the northern and western parts of the velocity map, labelled I1 to I5, named in clockwise order from the north.

Building on this segmentation, we defined a series of central positions along the isolated structures, spaced by ~100 pc. Around each centre, a localized VSF was computed from spaxels within a ~0.39 kpc radius aperture – chosen to match the widest filament widths in projection and to align with the MEGARA FOV. The results are not sensitive to the choice of radius, although changing it may shift the dominant turbulence scales being probed (see Sections 3.2.4 and 3.2.6). From each VSF, we extract three key parameters to quantify turbulence: (1) the slope (power-law index), measured from the seeing limit up to the flattening – or until velocity-pair statistics become insufficient if no flattening occurs; (2) the flattening scale, defined as the first point where the VSF flattens, which may indicate the characteristic injection scale of turbulent energy; and (3) the amplitude, measured as the VSF velocity value at the flattening, providing an estimate of largest velocity associated with the dominant turbulent motions. Flattening scale and amplitude are only measured when flattening is present. The same parameters were also extracted from MEGARA VSFs for the SE (1 comp.) and FE filaments to directly compare central filaments from MUSE with those at

larger projected distances. The full set of localized VSF regions and parameter variations – slope, flattening scale, and amplitude – as a function of distance from the galaxy centre (α: 12h30m49.4233s, δ: +12d23m28.043s; NED) and relative angle to the visible jet (oriented at 21.6°, measured from the projected jet axis in the radio image) are shown in Fig. 7. While substructures are defined to reduce projection contamination, they are not analysed individually; VSFs show broadly similar properties across those substructures, as seen in the following paragraph (with exception of F4).

The extracted VSF parameters exhibit spatial trends across the filament network. In particular, both slope and amplitude show a decline within the central ~0.8 kpc and remain generally lower at larger projected distances and angular separation from the jet axis. In the centre, slopes reach up to ~1.1 dex decreasing below 1 dex beyond the aperture diameter. Similarly, VSF amplitudes peak at ~100 km s⁻¹ near the nucleus and fall below ~50 km s⁻¹ further out. Two groups of localized outliers are highlighted in magenta: one projected near the jet axis (filled) and one associated with the subregion F4 (unfilled). Both exhibit elevated slope and amplitude values; the filled outlier matches central measurements, while the unfilled one is slightly higher. In contrast, the flattening scales show weaker trends with galactocentric distance or jet orientation. Most are below 0.2 kpc, though a few reach ~0.5 kpc, including the marked region near the jet axis. At larger distances, SE and FE filaments broadly follow the outer radial trend, with the FE filament showing slightly lower values. With respect to jet orientation, trends in VSF parameters are less pronounced. A modest decline from 0° to 90° may be present, but the highest values correspond to the innermost regions near the nucleus. The two outlier groups show distinct behaviours: the filled outlier close to the jet axis could reflect a genuine angular dependence, while the unfilled outlier deviates from other regions at similar position angles, likely due to local limitations (see Section 3.2.4). Overall, variations with angle appear to be shaped more by localized conditions than by a systematic geometric trend. Some southeastern filaments may also lie on the far side of the system, along the counter-jet direction not visible in observations (W. B. Sparks et al. 1993).

### 3.2.4 Interpretation of localized VSFs

Localized VSFs across the M87 filament network reveal both coherent trends and notable deviations. A particularly sharp decrease in VSF slope and amplitude occurs within ~0.8 kpc of the galaxy centre. These enhanced values are primarily associated with regions where the VSF aperture overlaps the nucleus, where complex velocity fields arise from the superposition of filamentary structures and the presence of a nuclear rotating ionized disc (see O + 23). In these cases, the resulting turbulence parameters are likely inflated due to superposition of distinct velocity components within a single aperture, rather than reflecting a purely physical change in the turbulence regime.

Beyond the central region – once the nuclear structure is excluded – the VSF parameters appear to stabilize. As seen from the median lines in the parameter plots, the slope, amplitude, and the flattening scale remain relatively consistent across most filamentary substructures, particularly in the central and SE filaments (single-component), and to a lesser degree in the FE filament. This uniformity suggests a common small-scale turbulence feature, possibly linked to a shared driver. A first candidate is energy injection from Type Ia supernovae (SNe Ia), which occur in the old stellar populations of massive ellipticals. Simulations (M. Li et al. 2020; R. Mohapatra & E.
Quataert (2024) show that SNe Ia drive compressive turbulence on ~100–200 pc scales with velocity dispersions up to ~40km s⁻¹ – consistent with the measured flattening scales and amplitudes. However, the greater distance of the FE filament from the main stellar body may call this explanation into question, though its slightly lower VSF parameters suggest they might not trace the same physical conditions. For this filament, this may reflect lower level of turbulence or structure formation, possibly linked to its location inside an older radio lobe associated with a past episode of AGN activity. In this scenario, the FE filament may be less actively driven, similar to the behavior on larger scales observed in outer filaments of many other galaxy clusters (Y. Li et al. 2020; S. Ganguly et al. 2023).

Alternative explanations for the observed flattening could also naturally account for the steeper slopes seen in several VSFs. Similar steep slopes, exceeding the classical Kolmogorov expectation, have been reported in other systems (e.g. Y. Li et al. 2020; S. Ganguly et al. 2023) and may arise when the VSF captures the influence of the primary driver (e.g. AGN feedback or SNe Ia) rather than the fully developed turbulent cascade. A possible explanation involves the partial dissipation of compressive modes on large scales prior to the transition to more purely Alfvénic turbulence on small scales with a Kolmogorov scaling, similar to what is observed in the solar wind plasma (H. Hu et al. 2022). Stratification may also affect the cascade across scales (R. Mohapatra et al. 2022; C. Wang, S. P. Oh & M. Ruszkowski 2023), although its influence at the scales probed here is likely modest. Other plausible contributors to the observed deviations include magnetic tension (C. J. Bambic, B. J. Morsony & C. S. Reynolds 2018; P. Grete, B. W. O’Shea & K. Beckwith 2021), plasma microinstabilities (M. W. Kunz, A. A. Schekochihin & J. M. Stone 2014; J. Squire, B. D. G. Chandran & R. Meyrand 2020; L. Arzamasskiy et al. 2023), gravitational acceleration (C. Wang et al. 2021), or ‘frozen-in’ supersonic turbulence (H. Hu et al. 2022). Altogether, these results suggest that local VSFs are highly sensitive to microphysical conditions and may reflect a superposition of regimes shaped by energy injection, dissipation processes, and magnetic fields.

Another factor affecting the interpretation of flattening features is the finite size of the analysis region. In both SE and FE filaments, flattening often appears below the half-height of the analysis window (aperture), but the measured flattening scale itself tends to correspond closely to half the width of the filament enclosed within the window. This suggests the flattening may partly reflect a geometric distortion related to the filament’s cross-section within the aperture, rather than a genuine physical dissipation scale. However, tests of the window size (see Section 3.2.6 and Appendix E) show that the flattening persists at similar scales even when the aperture is increased and include filaments of varying, mostly larger, widths.

Among the localized VSFs, a few outliers stand out with notably elevated slope and amplitude values (highlighted in magenta in Fig. 7). One of these lies near subregion F4, where the aperture covers overlapping filaments along the LOS – leading to artificially inflate turbulence measures due to projection confusion. VSFs marked with crosses in Fig. 7 (panel a), showing no flattening, are also likely affected by the blending of kinematically distinct filaments (e.g. F2 and F5 near F4) or by nuclear contamination in F3. The second outlier lies northeast of the nucleus, near the jet axis and within a filament bend. This region shows a high LOS velocity (~−275 km s⁻¹) and broad velocity dispersions (> 250 km s⁻¹ in the map from B + 19), in contrast to neighbouring areas, where velocities typically range from −50 to −100 km s⁻¹ and dispersion remain below 150 km s⁻¹. These features may indicate genuine jet interaction

or a geometrically thick, possibly LOS-elongated filament segment. Such a 3D extension would increase LOS velocity spread, thereby boosting the VSF amplitude relative to adjacent regions. In contrast, the SE filament – where distinct velocity components are more clearly resolved – shows a VSF flattening only when both components are blended in a single fit (1 comp.), whereas no flattening appears in the VSF when computed using only the first component. This suggests that some observed flattenings elsewhere may arise from unresolved LOS superposition rather than true physical injection or dissipation scales, underscoring the importance of disentangling kinematic substructure in multiphase turbulence analysis.

Finally, additional insight comes from comparing VSFs of the ionized and molecular gas in the SE filaments. Within the CO(2–1) detection region, both Hα- and CO-based VSFs exhibit a similar power-law slope at smaller scales. This behaviour aligns with multiphase turbulence models, where cold, dense gas inherits the motions (M. Gaspari & A. Sądowski 2017; M. Gaspari et al. 2018; G. M. Voit 2018; R. Mohapatra et al. 2022), and may also reflect direct coupling or dynamical interplay between phases at these scales. However, at larger separations (> ~0.15 kpc), the VSFs diverge: the CO VSF steepens then flattens near ~0.3 kpc, while Hα maintains a more consistent slope. Notably, when the Hα VSF is computed only within the CO detection footprint, it remains roughly parallel to the full SE-filament VSF but with a slightly lower amplitude. The divergence between the CO and Hα VSFs within the CO-detected region may arise from a combination of factors – including the differences in local volume-filling factor and projection effects. Additionally, the flattening of the CO VSF occurs near a scale comparable to half the width of the CO-emitting region (~300 pc east-west), suggesting that the finite spatial extent may also influence the VSF behaviour (see Section 3.2.6). Furthermore, to test whether the Hα within the CO detection VSF behaviour depends on aperture placement, we also computed it within a shifted CO-sized region elsewhere in the SE filaments. The small-scale slope is likely a robust feature, remaining consistent across apertures, while at larger scales (around the 1/2 window), the VSF steepens or flattens depending on location – suggesting greater sensitivity to geometry and local structure.

### 3.2.5 Turbulent heating in the FE filament

The dissipation of AGN-driven turbulence may serve as an important heating mechanism in cool-core clusters. In particular, I. Zhuravleva et al. (2014) showed that in M87 and Perseus, radiative cooling might be balanced by turbulent dissipation inferred from X-ray surface brightness fluctuations on large scales. However, X-ray observations lack the spatial resolution to probe turbulence on the smaller scales where most dissipation is expected to occur. S. Ganguly et al. (2023) extended this analysis to smaller scales by using Hα filaments as tracers of turbulence in the central regions of a few systems, including M87 with MUSE WFM data. Here, we extend their approach to larger radii by analysing the FE filament of M87 using MEGARA data.

We estimate the turbulent heating rate, $Q_{\text{turb}}$, at the seeing scale (1 arcsec ≈ 78 pc), following the formalism of S. Ganguly et al. (2023) for VSFs steeper than the Kolmogorov slope

$$Q_{\text{turb}} \approx 5.2 \frac{\rho v_{\text{l,1D}}^3}{l}, \qquad (6)$$

where $v_{\text{l,1D}}$ is the LOS velocity at scale $l$ from the VSF, and $\rho$ is the gas mass density given by $\rho = \mu m_p (n_i + n_e) = \mu m_p \xi n_e$. Here, $\mu = 0.61$ is the mean molecular weight, $m_p$ the proton mass, $\xi = 1.912$ accounts for total particle number (assuming 0.5 solar metallicity), and $n_e$ is the electron number density $n_e$. We compare this to the
radiative cooling rate

$$Q_{\text{cool}} = n_e^2(\xi - 1)\Lambda(T), \quad (7)$$

where $\Lambda(T)$ is the normalized cooling function computed from K. M. Schure et al. (2009), assuming solar metallicity. At ICM temperatures ($10^7 - 10^8$ K), the metallicity dependence is weak and can be reasonably neglected.

Adopting $n_e \approx 0.015 \text{ cm}^{-3}$ and $T \approx 1.9$ keV at 15 kpc (I. Zhuravleva et al. 2014), we estimate a turbulent heating rate of $Q_{\text{turb}} \approx 4.9 \times 10^{-28} \text{ erg s}^{-1} \text{ cm}^{-3}$ and a radiative cooling rate of $Q_{\text{cool}} \approx 4.1 \times 10^{-27} \text{ erg s}^{-1} \text{ cm}^{-3}$ for the FE filament. These results are consistent with findings from the central regions of M87 and other systems analysed by S. Ganguly et al. (2023), as well as with the behaviour observed in numerical simulations (e.g. C. S. Reynolds, S. A. Balbus & A. A. Schekochihin 2015; H. Y. K. Yang & C. S. Reynolds 2016; Y. Li, M. Ruszkowski & G. L. Bryan 2017), where turbulent heating estimated on small scales typically accounts for only a minor fraction of radiative cooling (here, ~10 per cent). This contrasts with earlier X-ray-based estimates (e.g. I. Zhuravleva et al. 2014, 2018) that probe larger scales and may overestimate the contribution of turbulent dissipation. A more detailed discussion of these differences and their implications can be found in S. Ganguly et al. (2023). We note that our estimated turbulent heating rate relies on the assumption that the kinematics of cool filaments and hot plasma are perfected coupled on small scales, which cannot be verified with existing X-ray telescopes yet due to their limited spatial resolution.

### 3.2.6 Limitations

It is important to note that the VSFs are subject to several limitations and sources of uncertainty that may affect the interpretation of results at various level. Below, we review key limitations based on conclusions drawn from previous studies and a few tests conducted on our data.

(i) **Projection effects.** Projection effects along the LOS introduce two primary biases in the derivation of VSFs. (1) LOS integration can blend uncorrelated velocity components, artificially lowering $\langle \Delta v \rangle$ at smaller scales. By contrast, side-by-side filaments retain their velocity contrasts in projection, boosting $\langle \Delta v \rangle$ at a given scale. Both effects are amplified in localized VSFs, where small apertures and limited spaxel coverage increase the likelihood of overlapping or adjacent structures dominating the measurement. This can steepen slopes and inflate amplitudes, particularly in regions like F4, where adjacent structures differ significantly in velocity (see Section 3.2.3). (2) Projection compresses true 3D separations into smaller 2D distances, shifting large velocity differences to artificially small scales and flattens the observed VSF. Numerical simulations (J. A. ZuHone, M. Markevitch & I. Zhuravleva 2016; R. Mohapatra et al. 2022; M. Fournier et al. 2025) demonstrate that projection steepens VSFs at small scales due to LOS cancellation but flattens them at larger scales where uncorrelated motions dominate. Although these competing effects cannot be disentangled without full 3D spatial information, some insight into LOS complexity can be gained from the spectral structure itself. In particular, multicomponent line profiles offer observational evidence for unresolved kinematic substructure along the LOS – as seen in other systems (e.g. N. A. Hatch et al. 2006; R. E. A. Canning et al. 2011; S. L. Hamer et al. 2016; B. Vigneron et al. 2024) – though their influence on VSFs has only recently been studied. C. Zhang et al. (2022) and M. Li et al. (2025) show that in compact sources, projection can flatten VSFs

by introducing high-|$\delta v|$ pixel pairs at small separations. Our results in the SE filaments (Section 3.2.1) support this: approximately 23 per cent of spaxels exhibit multiple velocity components, and the ~0.3 kpc flattening vanishes when the VSF is computed using only the first component. Moreover, the slope becomes slightly steeper in the single-component case, reinforcing that LOS superposition can mask true velocity correlations and bias turbulence estimates. We note that while MEGARA's resolution ($R \sim 6000$) allows us to resolve multiple components, finer substructure may remain unresolved. Higher resolution data could better constrain VSFs and reveal additional complexity.

(ii) **Seeing effects.** Atmospheric seeing impacts VSFs, particularly at small scales. Using MUSE data of M87, Y. Li et al. (2020) applied Gaussian smoothing matched to the point spread function (PSF) and found a modest (<0.1 dex) suppression of VSF power due to velocity homogenization between neighboring spaxels. Similar studies (M. C. Chen et al. 2023; S. Ganguly et al. 2023) report that smoothing steepens the VSF slopes below the PSF scale while preserving the main flattening. Our MEGARA-based tests using the same approach show similar suppression – despite fewer velocity pairs due to narrower filaments. Simulations by M. Fournier et al. (2025), modelling AGN-driven turbulence with the XMAGNET (eXascale simulations of Magnetized AGN feedback focusing on Energetics and Turbulence) framework, further confirm that smoothing artificially steepens VSF slopes below the kernel scale. Therefore, we consider the VSFs unreliable below the seeing limit.

(iii) **Spaxel size and binning.** Our data sets with span different spaxel sizes, from 0.2 arcsec (MUSE) and 0.31 arcsec (MEGARA) to 0.64 arcsec for SITELLE (after $2 \times 2$ binning), raising concerns about their effect on VSFs. To evaluate this, we conducted a series of spaxel size/binning tests using SITELLE, MEGARA, and MUSE data. As detailed in Appendix D, we compared matched subregions using the original SITELLE spaxel size and both original and resampled MUSE data. We also applied various binning schemes to MUSE (up to $10 \times 10$) and MEGARA ($3 \times 3$). In all cases, the overall VSF shapes remain relatively stable. Larger spaxel sizes occasionally steepen the slope at small scales, but flattening scale and amplitude stay consistent. These results align with simulations (e.g. M. Fournier et al. 2025), which report that resolution degradation can bias small-scale VSF properties. Although such effects are relatively minor within our tested spaxel range, they should none the less be kept in mind when interpreting VSFs. Notably, SITELLE consistently shows higher amplitudes than MUSE or MEGARA, even after matching spaxel size. We therefore attribute this difference to the noise of SITELLE's velocity map, which introduces larger spaxel-to-spaxel variations and enhances apparent velocity differences.

(iv) **Window size.** To evaluate whether the size and geometry of the analysis window affects VSF measurements, we used the full MUSE velocity map of the central filaments and computed VSFs within circular and rectangular windows of varying sizes. While only the circular-window results are shown in Appendix E, we find that both geometries yield consistent VSF shapes. Across window sizes, the first flattening scale and amplitude remain stable, suggesting that small- and intermediate-scale dynamics are unaffected. Differences mainly appear at larger separations, where a second flattening or a diminished first flattening can emerge due to contributions from large-scale dynamics. This implies that window size does not strongly bias VSFs within the explored FOV range. However, larger windows naturally incorporate broader scales, potentially shifting the dominant dynamics from local turbulence to global kinematics. Few studies (e.g. T. Ha et al. 2022; M. Li et al. 2025) have tested the impact of window size on VSFs, showing that beyond separations of roughly
half the map size, edge effects, and reduced pair statistics may lead to artificial flattenings. Although our own test yields broadly consistent trends, this underscores the need for caution when interpreting VSFs at the largest scales.

(v) Noise. The impact of noise on VSFs has been explored in previous work. M. Li et al. (2020) demonstrated that higher noise levels, modelled through lower SNR cuts, can raise the VSF amplitude at small scales and slightly alter the slope, while leaving the flattening unaffected. This may partly explain the systematically higher VSF amplitudes seen in SITELLE, likely driven by increased noise and its lower spectral resolution.

# 4 IONIZATION SOURCES

## 4.1 Line ratios

Several emission lines in the MEGARA and SITELLE data – specifically [O I]λ6300, Hα, [N II]λ6583, and [S II]λ6716, 6731 – provide diagnostics of the physical conditions of the ionized gas. As a first step, we examine key emission-line ratios sensitive to ionization mechanisms such as shocks: [O I]/Hα$^6$, [N II]/Hα, and [S II]/Hα.$^7$

Due to the faintness and low SNR of the [O I]λ6300 line in MEGARA spaxels, spatially resolved [O I]/Hα maps could not be constructed. Instead, we measure the [O I]/Hα ratio from the integrated MEGARA spectra of the SE and FE filaments, obtaining low values of 0.23 and 0.16, respectively. These results are consistent with previous measurements based on MUSE WFM data for the SE filaments (~0.2; B + 19) and long-slit spectroscopy from the Observatoire de Haute Provence for the FE filament (0.18; G + 00).

The [N II]/Hα and [S II]/Hα ratios are sufficiently strong to be mapped with SITELLE and MEGARA. Fig. 8 presents these ratio maps along with radial gradients derived from SITELLE data, computed using concentric 1 kpc-wide annuli centred on the galaxy nucleus. Tests confirmed that the analysis is insensitive to annulus spacing. Within each annulus, line ratio distributions are summarized with boxplots (median, interquartile range, and overall spread), and a weighted linear regression on median values to quantify the radial trends. Overall, SITELLE results are consistent with MEGARA and MUSE WFM but systematically slightly higher, likely due to increased noise levels. In the central regions, SITELLE shows high [N II]/Hα ratios of ~2. Beyond ~2.5 kpc, the gradient does not decline monotonically: local variations are observed between 2.5–5 kpc, especially toward the southern filament. Farther out, the ratio gradually decreases to ~1.5 at ~7.5 kpc. A similar, though shallower, trend is seen in the [S II]/Hα profile, decreasing from ~1.3 to ~1.1, with greater dispersion at intermediate radii. MEGARA observations of the SE filaments confirm the elevated ratios seen with SITELLE. The [N II]/Hα ratio decreases from ~2.75 near the nucleus to ~1.95 further out, with the highest values projected along the edge of the radio lobe. This regions also corresponds to where two velocity components are detected, both exhibiting similar elevated values around ~2.5. The [S II]/Hα ratio follows a similar pattern, declining from ~1.51 to ~1.24, with both components near ~1.45. N. Werner et al. (2013) and B + 19 similarly reported high [O III]/Hβ ratios in this same region (~2) – exceeding those measured in the galaxy nucleus – consistent with Seyfert-like ionization. By contrast, the FE filament displays significantly lower ratios in both data sets. The [N II]/Hα ratio remains near ~1.0, and the [S II]/Hα ratio near ~0.5,

consistent with a lower ionization environment. Although we do not measure the [O III] and Hβ lines in our new data, we recall that G. Gavazzi et al. (2000) reported a [O III]/Hβ ratio of ~0.85 in the same region.

[S II]λ6716/[S II]λ6731 ratio maps were also produced from the MEGARA data (see Fig. 8), serving as a tracer of the electron density $n_e$ of the ionized gas (D. E. Osterbrock & G. J. Ferland 2006). This analysis was not performed with the SITELLE data due to its lower spectral resolution and insufficient SNR for [S II]λ6731 across the field. In the SE filaments, the mean [S II] doublet ratio is ~1.37, corresponding to $n_e$ ~ 40 cm⁻³ assuming an electron temperature of 10 000 K and using the calibration from B. Proxauf, S. Örtl & S. Kimeswenger (2014). This ratio is consistent with previous measurements, including ~1.40 reported by N. Werner et al. (2013) and B + 19 at the edges of the SE filaments. In the FE filament, we measure a higher mean ratio of ~1.56, which lies above the low density limit with $n_e$ ~ 10 cm⁻³, as reported by G + 00.

## 4.2 WHAN diagram

To further investigate the nature of the ionizing sources in the SE and FE filaments, we use the width Hα-nitrogen (WHAN) diagram (R. Cid Fernandes et al. 2010, 2011), traditionally used to classify galaxies by dominant ionization mechanisms. It an alternative to the classical Baldwin, Phillips & Terlevich diagram (J. A. Baldwin, M. M. Phillips & R. Terlevich 1981) when key emission lines such as [O III]λ5007 and Hβ are unavailable, as in our MEGARA and SITELLE data. However, these lines have been detected with other instruments in the central (B + 19), southern (N. Werner et al. 2013) and in the outer FE filament (G + 00). The WHAN diagram combines the [N II]/Hα line ratio with the Hα equivalent width $W_α$, defined as line flux over local continuum flux, measuring line strength relative to the stellar background.

In this framework, we use well-established diagnostic lines to guide our interpretation. Specifically, regions with $log([N II]/Hα) > 0.4$ and $W_α > 6$ Å are classified as Seyfert-like, indicative of strong ionizing sources such as AGN; those with $3$ Å < $W_α$ < 6 Å fall into the LINER (low-ionization nuclear emission-line region)-like regime, typically associated with weaker ionization. Values of $W_α < 3$ Å correspond to retired or passive regions, where ionization is likely driven by old stellar populations. Star-forming regions are found at $log([N II]/Hα) < -0.4$ and $W_α > 3$ Å, consistent with photoionization by young massive stars. Here, we use the WHAN diagram for the SE and FE filaments – not to impose a strict classification, but to identify regions of relatively higher or lower energy ionization source, acknowledging the likely contribution of multiple excitation mechanisms (G. J. Ferland et al. 2009; M. McDonald, S. Veilleux & D. S. N. Rupke 2012; C. L. Rhea et al. 2025).

Fig. 9 presents the WHAN diagram on a spaxel-by-spaxel basis for the SE and FE filaments, along with the equivalent width maps. In the SE region, where two components were identified, both components are shown separately. No spaxels fall in the star-forming WHAN region, consistent with the high [N II]/Hα ratios presented in Section 4.1. The FE filament primarily occupies the Seyfert-like region, with a mean $W_α$ ~ 150 Å. In contrast, the SE filaments initially show much lower $W_α$ values (mean ~5 Å) due to dilution by the underlying stellar continuum, placing spaxels in the passive, LINER-like, or Seyfert-like regions.

We note that both filaments exhibit broad spread in $W_α$. Measurement uncertainties may contribute but geometry– including projection effects, varying thickness, partial LOS overlap – likely

<footnote>$^6$[O I]λ6300.</footnote>

<footnote>$^7$([S II]λ6716 + [S II]λ6731)/Hα.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Maps of the [N II]/Hα (top panel), [S II]/Hα (middle panel), and [S II]λ6716/[S II]λ6731 (bottom panel) line ratios, using only spaxels with SNR<sub>Hα</sub> > 3 and SNR<sub>[SII]λ6716</sub> > 3 where applicable, derived from (a) MEGARA and (b) SITELLE observations. In the MEGARA maps, black gradient contours over the SE filament indicate the 6 cm VLA radio emission tracing the counter-jet lobe at 3σ and 10σ (D. C. Hines et al. 1989). In the SITELLE maps, the MEGARA pointings targeting the SE and FE filaments are outlined in white boxes. (c) A radial gradient was also derived for the two first line ratios with SITELLE data, using concentric annuli centred on the galaxy nucleus. Distributions of the line ratio within each annulus are shown as boxplots, from which the median values were used to perform a linear regression. The best-fitting parameters of the regression are reported in the upper right corner of the panel. Colour scales are consistent across maps of the same type. All scale bars correspond to 5 arcsec.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Maps of the $H\alpha$ equivalent width $W_\alpha$ for the SE and FE filaments observed with MEGARA ($\text{SNR}_{H\alpha} > 3$). The white outline marks the two-component region within SE filaments. (b) WHAN diagrams constructed spaxel-by-spaxel with each filament (SE in pink – 1 comp., right; and FE in orange, left). The mean values for each distribution are marked by a black square (SE) and a black circle (FE). (c) WHAN diagram showing the first and second components within the SE two-component region. Both components largely overlap with the single-component SE distribution (black contour).</figcaption></figure>

plays a role. To assess measurement impact, we examined $W_\alpha$ versus $H\alpha$ SNR, finding strong positive correlations (Pearson coefficients 0.93 for SE, and 0.76 for FE). This sensitivity highlights limitations of spaxel-level WHAN diagnostics and suggests integrated-region analyses yield more reliable ionization classifications.

## 4.3 Velocity dispersion versus line ratios

Combining ionization-sensitive line ratios with kinematics can also offer insight into the excitation mechanisms. In particular, diagrams comparing line ratios to velocity dispersion are useful to identify shocks, which tend to broaden lines and enhance low-ionization transitions such as [N II], [S II], and [O I] relative to $H\alpha$ (J. A. Rich, L. J. Kewley & M. A. Dopita 2011, 2015; M. McDonald et al. 2012; I. T. Ho et al. 2014; A. Boselli et al. 2019; V. Johnston 2022). In contrast, photoionized regions typically exhibit both lower line ratios and narrower velocity dispersions.

In M87, B + 19 identified a clear positive correlation between velocity dispersion and [O I]/$H\alpha$ in the central filaments, consistent with AGN-driven shocks, with weaker or absent correlations for [N II]/$H\alpha$ and [S II]/$H\alpha$ at larger distances. Here, we extend this analysis to the SE and FE filaments using MEGARA data to assess whether similar trends persist farther from the nucleus. We constructed spaxel-by-spaxel diagrams of velocity dispersion versus [N II]/$H\alpha$ and [S II]/$H\alpha$ for each pointing and velocity component (where applicable), and quantified the trends using the Pearson correlation coefficient *r*. These results are presented in Fig. 10.

No strong or significant correlations emerge in either filament. [N II]/$H\alpha$ generally yielded slightly higher *r* values than [S II]/$H\alpha$. Overall, both filaments show relatively high line ratios but low

velocity dispersions – typically <100 km s⁻¹, and only a few tens of km s⁻¹ in the FE filament, consistent with earlier results (G + 00).

## 4.4 Interpreting the ionization mechanisms

The spatial and spectral properties of warm ionized gas in M87 reveal a complex excitation environment with multiple ionization processes contributing. Our results support previous findings that photoionization by young massive stars is unlikely to dominate, given the absence of compact UV or optical stellar counterparts associated with the filaments (N. Werner et al. 2013, B + 19; P. Tamhane et al. 2025) and the elevated low-ionization line ratios observed throughout the system. The decoupling between the ionized gas and the stellar continuum, particularly in the FE filament, further rules out stellar sources and instead implies the presence of a harder excitation mechanism. Although, we note that post-asymptotic giant branch stellar populations can produce LINER-like emission (e.g. P. Lagos et al. 2022). Such line ratios are primarily observed in the central region (as shown in B + 19) and in the SE filaments.

Ionization conditions vary spatially. Within the central ~0.5 kpc, [N II]/$H\alpha$ and [S II]/$H\alpha$ ratios are high, coinciding with broad velocity dispersions and previously reported elevated [O I]/$H\alpha$ and [O III]/$H\beta$ line ratios from MUSE (B + 19). Collectively, these diagnostics point toward excitation dominated by the AGN.

The SE filaments exhibit similarly high, or even higher, line ratios – especially near the projected boundary of the counter-jet radio lobe. Ratios decline sharply beyond the region where two velocity components are detected, just outside the radio lobe' edge.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Spaxel-by-spaxel diagrams comparing the velocity dispersion to (a) [N II]/Hα and (b) [S II]/Hα for the SE and FE filaments from MEGARA data. Data points from the SE filaments are shown in light grey circles (1 comp.), dark grey squares (1st comp., in the two-component region), and pink triangles (2nd comp.). The FE filament is represented by black hexagons. The Pearson correlation coefficient *r* is indicated in the upper left corner or each panel; values near ±1 indicate strong positive/negative correlation, while values near 0 indicate none.</figcaption></figure>

This spatial correspondence suggests direct mechanical interaction between the AGN outflow and the filamentary gas, likely via shocks or compressive flows. However, we find no clear correlation between line ratios and velocity dispersion – a classical signature of fast radiative shocks. Similar results have been reported in other BCGs (e.g. S. L. Hamer et al. 2016, 2019), where high-excitation line ratios lack the broad profiles predicted by shock models. Moreover, the observed line ratios – particularly [O I]/Hα and [S II]/Hα – disagree with the predictions of fast shock models (e.g. M. G. Allen et al. 2008). The relatively low velocity dispersions measured here suggest either weak shocks or additional ionization sources. Plausible alternatives include collisional heating by cosmic rays, potentially energized by magnetohydrodynamic waves along the filaments (G. J. Ferland et al. 2009; A. C. Fabian et al. 2011; M. McDonald et al. 2012), dust scattering (K. Mattila, M. Juvela & K. Lehtinen 2007; K.-I. Seon & A. N. Witt 2015), thermal conduction from the surrounding hot ICM, and ionization by soft X-ray/extreme UV (EUV) photons from the cooling ICM (F. L. Polles et al. 2021).

In contrast, the FE filament has a lower ionization state, with significantly lower line ratios and narrow line widths. None the less, the elevated Hα equivalent widths, persistent [N II]/Hα ratio, and lack of underlying stellar continuum suggest that stellar photoionization is not the dominant excitation mechanism. G + O0 interpreted this outer filament as relic of past AGN activity, potentially energized by intermediate-velocity shocks from the now-inactive eastern radio lobe – a scenario consistent with its morphology, location, and ionization properties. Additionally, the [O I]/Hα ratio again disagrees with fast shock models (M. G. Allen et al. 2008), and the absence of broadened lines may point to fossil or decaying shocks, or steady-state processes such as cosmic ray heating, turbulent mixing layers operating at lower energies, or X-ray emitting gas from the cooling ICM (e.g. A. C. Fabian et al. 2011; M. McDonald et al. 2012; S. L. Hamer et al. 2016; F. L. Polles et al. 2021).

## 5 Hα-X-RAY CORRELATION

A spatial correlation between Hα + [N II] emission and soft X-ray structures in M87 has long been noted, particularly within the central ~5 kpc (A. J. Young et al. 2002; W. B. Sparks et al. 2004; N. Werner et al. 2010, 2013). These studies highlighted a close morphological association between warm ionized gas and the cooler X-ray-emitting phase (~0.5 keV), suggesting thermal coupling and/or co-evolution. More recently, V. Olivares et al. (2025) quantified a strong correlation between Hα and X-ray surface brightness in seven cool-core clusters, including M87 using the MUSE WFM data covering the central ~1 arcmin × 1 arcmin region. In this work, we revisit this correlation using SITELLE's larger FOV, covering the southeastern and FE filaments beyond the central region.

We adopt the same methodology as presented in V. Olivares et al. (2025) to compare the Hα and soft X-ray surface brightness in M87. We use deep archival *Chandra* observations of M87, from which the filamentary X-ray component is isolated using the Poisson Generalized Morphological Component Analysis (pGMCA; J. Bobin et al. 2020). The pGMCA method can extract a user-defined number of components with no prior physical or instrumental knowledge, by focusing on the morphological diversity of the different emissions. We decompose the X-ray emission into four components reflecting a mix of structures: diffuse halo, AGN cavities, and filaments – without a strict one-to-one association. In our case, the filamentary emission is captured across two distinct components. The first highlights the sharper, well-defined filaments, broadly consistent with the morphology of the warm ionized gas. The second component also traces some of these filamentary structures, but additionally includes more extended regions lacking a clear Hα counterparts – notably encompassing parts of the X-ray arms (e.g. southwestern arm). For the warm ionized phase, we use SITELLE Hα map (binned 2 × 2). Filamentary regions are identified manually on this map, following
<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Hα flux map (SNR$Hα > 3, 2 \times 2$ binned). (b) X-ray flux map of the first component from the decomposition. (c) X-ray flux map of the second component. In panels (a)–(c), surface brightness extraction regions are shown in white/red (same regions across panels), and the MUSE WFM pointing is indicated by the light grey dotted box. (d) Hα–X ray surface brightness correlation for filaments in M87 (red circles, representing the sum of both X-ray components) and two isolated clumps in Centaurus (blue squares). For M87, crosses and white circles with red contours represent values derived using the first and second X-ray components, respectively. Errors are smaller than the scatter. The dashed black line shows the best-fitting relation from V. Olivares et al. (2025), and the thick grey line corresponds to the diffuse gas trend in stripped tails from M. Sun et al. (2021). (e) Full FOV X-ray map of the second component, including the southwestern arm. The zoomed-in region shown in panels (b) and (c) is indicated by the box.</figcaption></figure>

coherent morphological structures in projection, and applied to the X-ray filament maps to ensure spatial correspondence. For each region, we extract both the Hα and soft X-ray fluxes, compute the luminosities using the redshift-dependent luminosity distance, and derive the surface brightness (in erg s⁻¹ kpc⁻²) by dividing by the physical area. Regions contaminated by the central AGN or background sources are excluded. Details of the surface brightness fitting and uncertainty treatment is provided in V. Olivares et al. (2025).

Fig. 11 shows the Hα–X-ray surface brightness correlation across M87’s filaments, using the two X-ray components identified as containing filamentary emission. The first and second components yield individual X-ray/Hα surface brightness ratios of 1.3<sup>+0.6</sup><sub>-1.1</sub> and 3.7<sup>+2.7</sup><sub>-2.3</sub>, respectively. Summed together, they reproduce the tight

correlation previously observed with MUSE WFM (V. Olivares et al. 2025), now extend to the larger filament network captured by SITELLE. The typical X-ray/Hα surface brightness ratio for the combined components (excluding the FE filament) is 5.0<sup>+3.4</sup><sub>-3.0</sub>, consistent with the cool-core average of 4.1±2.4 found in that study. This consistency reinforces the idea that local excitation mechanisms (e.g. cosmic ray heating, turbulent mixing layers, and reprocessed X-ray radiation) govern the emission, while the distance to the AGN plays little role. The persistence of a strong correlation across diverse environments, including stripped galaxy tails (M. Sun et al. 2021), supports a common mechanism linking hot and warm gas. While this behaviour is consistent with models in which local thermal instabilities in the ICM lead to condensation of warm gas (M. McCourt et al. 2012; P. Sharma et al. 2012; Y. Li & G. L. Bryan 2014),
it may more broadly reflect ongoing dynamical interactions between warm filaments and the surrounding hot atmosphere – regardless of whether the warm gas formed via precipitation or uplift-driven processes.

Intriguingly, one region in M87—the isolated FE filament (labelled ‘FE’ in Fig. 11)—lies well below the main correlation with X-ray/Hα ratio of 0.35. A similar deviation is found in Centaurus cluster, where we applied the same pGMCA-based analysis with four X-ray components. While V. Olivares et al. (2025) found that Centaurus’ central filaments follow the main trend, our extended analysis identifies two more isolated northeastern clumps. Among these, the clump located at 7–9 kpc (clump 1 in panel d of Fig. 11) presents a low X-ray/Hα ratio of 1.2. In both Centaurus and M87, the deviating clumps lie near the inner edges of extended soft X-ray features likely shaped by past AGN activity: the northeastern Centaurus clump aligns with the direction of a ghost cavity (C. S. Crawford et al. 2005; J. S. Sanders et al. 2016), while M87’s FE filament sits within an old radio lobe. This suggests that such clumps were uplifted during earlier episodes of radio-mode feedback and have since cooled more, evolved differently or become less coupled to the surrounding hot ICM compared to the core-connected filaments. Alternatively, the observed discrepancy may reflect limitations in isolating diffuse filamentary X-ray emission using pGMCA: as the algorithm can extract components based on their morphological features, it will more easily retrieve sources that are extended with a characteristic shape than small clumps – particularly in low surface brightness regions.

## 6 SUMMARY AND CONCLUSIONS

We presented a comprehensive analysis of the warm ionized filaments in M87, combining new and archival multiwavelength data to probe their kinematics, turbulence, ionization state, and multiphase connections in different environments shaped by AGN feedback. MEGARA observations target two key regions: the SE filaments, which host the only known molecular cloud detected with ALMA, and the isolated FE filament within an old radio lobe. SITELLE offers full coverage of the filamentary system, enabling a global view of the ionized gas. These data sets were complemented by MUSE WFM observations of the central filaments and *Chandra* X-ray data tracing the surrounding hot plasma. Our main findings are as follows:

(i) M87’s ionized filaments exhibit complex motions that reflect strong local influences (e.g. jet interaction) beyond simple rotation. SITELLE maps the entire velocity field for the first time, reaching outer filaments skirting buoyant bubbles of different ages with velocities from ~−250 to 50 km s⁻¹. MEGARA resolves a two-component structure in the SE filaments near the radio lobe edge, separated by ~70 km s⁻¹, with similar velocity dispersions of ~30–50 km s⁻¹. The FE filament shows lower velocities and narrower dispersion (~60 km s⁻¹, and σ ~30 km s⁻¹).

(ii) We computed VSFs across localized subregions of the filament network. Most VSFs show steep slopes (> ~2/3) and flattening at small scales of a few hundred parsecs. These trends, consistent across the central and SE filaments, point to a common turbulence driver – most likely SNe Ia and, if confirmed, observed here for the first time – rather than the cascade of turbulence and/or shaped by dissipation processes or magnetic effects. In contrast, the FE filament exhibits slightly lower amplitude, suggesting weaker active driving. The proximity of M87 allows us to resolve these scales, but higher resolution observation (e.g. MUSE NFM) are needed to probe smaller scales and to assess the robustness of these findings.

In some regions, elevated VSF amplitudes or slopes are observed, likely due to LOS superposition. We also emphasize that the finite size of the analysis apertures can influence the measured flattening. Such limitations must be considered when interpreting the observed turbulence cascade.

(iii) In the SE filament, CO(2–1) emission detected with ALMA coincides with the ionized gas, and both phases show comparable VSFs on small scales. This agreement is consistent with models where colder gas inherits warm-gas motions, and may reflect direct dynamical coupling between the phases. Differences at larger scales and broader ionized-gas velocity dispersions may reflect geometry and differing volume filling factors.

(iv) Emission-line diagnostics, including line ratios and WHAN diagrams, indicate that ionization in the filaments is dominated by non-stellar sources. The elevated [N II]/Hα and [S II]/Hα ratios, especially near the radio lobes, are consistent with AGN-related processes such as slow shocks, and may also arise from cosmic ray heating, turbulent mixing layers or X-ray emitting gas from the cooling ICM that is triggered by AGN feedback. The FE filament shows lower excitation but exceeds levels expected from old stellar populations, consistent with fossil AGN activity or low-energy processes.

(v) We confirm a strong spatial correlation between Hα and soft X-ray surface brightness in the main filaments, as seen in other cool core clusters, implying close coupling between warm and hot phases. This behaviour may result from thermal instabilities in the ICM or more generally from the interaction between cold gas and the surrounding hot atmosphere. However, this correlation weakens in the isolated FE filament and in an outer clump in Centaurus, both with lower X-ray/Hα ratios (~1), possibly representing to gas uplifted during earlier AGN activity that has since cooled more efficiently, evolved differently, or become less coupled to the ambient hot phase. However, we note that this discrepancy could also reflect pGMCA limitations.

Overall, this study demonstrates the value of combining spatially resolved spectroscopy with multiphase to investigate AGN-driven filamentary structures in the core of M87. By targeting both central and more distant regions (e.g. FE filament), we highlight how local conditions and feedback history shape the dynamics, excitation, and coupling of the multiphase gas.

## ACKNOWLEDGEMENTS

We thank GTC staff for their flexibility (allowing visitor mode observations), help, and support for data acquisition, quality assessment, and reduction (D. Reverte Payá), as well as co-observer R. Varas González for his support during the GTC observing run. The authors acknowledges that this work is based on observations obtained with SITELLE, a joint project of Université Laval, ABB, Université de Montréal, and the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. The observations were carried out with care and respect from the summit of Maunakea, a site of great cultural and historical significance. CP acknowledges support from the Fonds de recherche du Québec-Nature et technologies (FRQNT) master's training scholarship (https://doi.org/10.69777/346450). M-LG-M acknowledges financial support from the grant CEX2021-001131-S funded by MCIU/AEI/ 10.13039/501100011033, from the coordination of the participation in SKA-SPAIN, funded by the
Ministry of Science, Innovation and Universities (MCIU), as well as NSERC via the Discovery grant program and the Canada Research Chair program. VO acknowledges support from the DYCIT ESO-Chile Comite Mixto PS 1757, Fondecyt Regular 1251702. YL acknowledges support from NASA grant 80NSSC22K0668, *Chandra X-ray Observatory* grant TM3-24005X, NSF grants AST-2514692 and AST-2510198, and CAREER award AST-2516092. LHM acknowledges finantial support by the grant PID2021-124665NB-I00 funded by the Spanish Ministry of Science and Innovation and the State Agency of Research MCIN/AEI/10.13039/501100011033 PID2021-124665NB-I00 and ERDF A way of making Europe. SC acknowledges the I+D + i project PID2022-140871NB-C21, financed by MICIU/AEI/10.13039/501100011033/ and 'FEDER/UE'.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the authors.

REFERENCES

Abdurro'uf et al., 2022, ApJS, 259, 35

Akaike H., 2011, in Int. Encyc. Stat. Sci. Springer, Berlin, Heidelberg, p. 25

Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J., 2008, ApJS, 178, 20

Arp H. C., 1967, Astrophys. Lett., 1, 1

Arroyo-Polonio A. et al., 2024, A&A, 687, A77

Arzamasskiy L., Kunz M. W., Squire J., Quataert E., Schekochihin A. A., 2023, Phys. Rev. X, 13, 021014

Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5

Bambic C. J., Morsony B. J., Reynolds C. S., 2018, ApJ, 857, 84

Blakeslee J. P. et al., 2009, ApJ, 694, 556

Bobin J., El Hamzaoui I., Picquenot A., Acero F., 2020, IEEE Trans. Image Process., 29, 9429

Boldyrev S. A., 1998, Phys. Plasmas, 5, 1681

Boselli A. et al., 2018, A&A, 614, A56

Boselli A. et al., 2019, A&A, 623, A52

Canning R. E. A., Fabian A. C., Johnstone R. M., Sanders J. S., Crawford C. S., Ferland G. J., Hatch N. A., 2011, MNRAS, 417, 3080

Cantiello M. et al., 2024, ApJ, 966, 145

Cappellari M., Emsellem E., 2004, PASP, 116, 138

Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245

Cazzoli S. et al., 2018, MNRAS, 480, 1106

Chen M. C. et al., 2023, MNRAS, 518, 2354

Churazov E., Brüggen M., Kaiser C. R., Böhringer H., Forman W., 2001, ApJ, 554, 261

Churazov E., Forman W., Vikhlinin A., Tremaine S., Gerhard O., Jones C., 2008, MNRAS, 388, 1062

Cid Fernandes R., Stasińska G., Schlickmann M. S., Mateus A., Vale Asari N., Schoenell W., Sodré L., 2010, MNRAS, 403, 1036

Cid Fernandes R., Stasińska G., Mateus A., Vale Asari N., 2011, MNRAS, 413, 1687

Crawford C. S., Hatch N. A., Fabian A. C., Sanders J. S., 2005, MNRAS, 363, 216

Drissen L. et al., 2019, MNRAS, 485, 3930

Emsellem E., Krajnovic D., Sarzi M., 2014, MNRAS, 445, L79

Event Horizon Telescope Collaboration, 2019, ApJ, 875, L1

Fabian A. C., Sanders J. S., Williams R. J. R., Lazarian A., Ferland G. J., Johnstone R. M., 2011, MNRAS, 417, 172

Fabian A. C., Ferland G. J., Sanders J. S., McNamara B. R., Pinto C., Walker S. A., 2022, MNRAS, 515, 3336

Falcón-Barroso J., Sánchez-Blázquez P., Vazdekis A., Ricciardelli E., Cardiel N., Cenarro A. J., Gorgas J., Peletier R. F., 2011, A&A, 532, A95

Ferland G. J., Fabian A. C., Hatch N. A., Johnstone R. M., Porter R. L., van Hoof P. A. M., Williams R. J. R., 2009, MNRAS, 392, 1475

Ford H. C., Butcher H., 1979, ApJS, 41, 147

Forman W. et al., 2005, ApJ, 635, 894

Forman W. et al., 2007, ApJ, 665, 1057

Forman W., Churazov E., Jones C., Heinz S., Kraft R., Vikhlinin A., 2017, ApJ, 844, 122

Fournier M., Grete P., Brüggen M., O'Shea B. W., Prasad D., Wibking B. D., Glines F. W., Mohapatra R., 2025, A&A, 698, A121

Gaia Collaboration, 2023, A&A, 674, A1

Ganguly S. et al., 2023, Front. Astron. Space Sci., 10, 1138613

Gaspari M., Sądowski A., 2017, ApJ, 837, 149

Gaspari M., Ruszkowski M., Oh S. P., 2013, MNRAS, 432, 3401

Gaspari M., Temi P., Brighenti F., 2017, MNRAS, 466, 677

Gaspari M. et al., 2018, ApJ, 854, 167

Gavazzi G., Boselli A., Vílchez J. M., Iglesias-Paramo J., Bonfanti C., 2000, A&A, 361, 1

Gebhardt K., Thomas J., 2009, ApJ, 700, 1690

Grete P., O'Shea B. W., Beckwith K., 2021, ApJ, 909, 148

Ha T., Li Y., Kounkel M., Xu S., Li H., Zheng Y., 2022, ApJ, 934, 7

Hamer S. L. et al., 2016, MNRAS, 460, 1758

Hamer S. L. et al., 2019, MNRAS, 483, 4984

Hatch N. A., Crawford C. S., Johnstone R. M., Fabian A. C., 2006, MNRAS, 367, 433

Heckman T. M., Baum S. A., van Breugel W. J. M., McCarthy P., 1989, ApJ, 338, 48

Hines D. C., Owen F. N., Eilek J. A., 1989, ApJ, 347, 713

Hlavacek-Larrondo J., Li Y., Churazov E., 2022, in Bambi C., Sangangelo A. eds, Handbook of X-ray and Gamma-ray Astrophysics. Springer, Singapore, p. 5

Ho I. T. et al., 2014, MNRAS, 444, 3894

Hu H., Qiu Y., Gendron-Marsolais M.-L., Bogdanović T., Hlavacek-Larrondo J., Ho L. C., Inayoshi K., McNamara B. R., 2022, ApJ, 929, L30

Johnston V., 2022, in American Astronomical Society Meeting #240, p. 230.01

Kolmogorov A., 1941, Akad. Nauk SSSR Dokl., 30, 301

Kunz M. W., Schekochihin A. A., Stone J. M., 2014, Phys. Rev. Lett., 112, 205003

Lagos P. et al., 2022, MNRAS, 516, 5487

Li M., Li Y., Bryan G. L., Ostriker E. C., Quataert E., 2020, ApJ, 898, 23

Li M. et al., 2025, ApJ, 984, 22

Li Y., Bryan G. L., 2014, ApJ, 789, 153

Li Y., Bryan G. L., Ruszkowski M., Voit G. M., O'Shea B. W., Donahue M., 2015, ApJ, 811, 73

Li Y., Ruszkowski M., Bryan G. L., 2017, ApJ, 847, 106

Li Y. et al., 2020, ApJ, 889, L1

Liepold E. R., Ma C.-P., Walsh J. L., 2023, ApJ, 945, L35

Macchetto F., Pastoriza M., Caon N., Sparks W. B., Giavalisco M., Bender R., Capaccioli M., 1996, A&AS, 120, 463

Martin T. B., Drissen L., Melchior A.-L., 2018, MNRAS, 473, 4130

Martin T., Drissen L., Prunet S., 2021, MNRAS, 505, 5514

Mattila K., Juvela M., Lehtinen K., 2007, ApJ, 654, L131

Mayer L., Mastropietro C., Wadsley J., Stadel J., Moore B., 2006, MNRAS, 369, 1021

McCourt M., Sharma P., Quataert E., Parrish I. J., 2012, MNRAS, 419, 3319

Donald M., Veilleux S., Rupke D. S. N., 2012, ApJ, 746, 153

McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117

McNamara B. R., Russell H. R., Nulsen P. E. J., Hogan M. T., Fabian A. C., Pulido F., Edge A. C., 2016, ApJ, 830, 79

Million E. T., Werner N., Simionescu A., Allen S. W., Nulsen P. E. J., Fabian A. C., Böhringer H., Sanders J. S., 2010, MNRAS, 407, 2046

Mohapatra R., Quataert E., 2024, ApJ, 965, 105

Mohapatra R., Jetti M., Sharma P., Federrath C., 2022, MNRAS, 510, 2327

Newville M., Stensitzki T., Allen D. B., Rawlik M., Ingargiola A., Nelson A., 2016, Astrophysics Source Code Library, record ascl:1606.014

Olivares V. et al., 2019, A&A, 631, A22

Olivares V., Su Y., Nulsen P., Kraft R., Somboonpanyakul T., Andrade-Santos F., Jones C., Forman W., 2022a, MNRAS, 516, L101

Olivares V. et al., 2022b, A&A, 666, A94

Olivares V. et al., 2023, ApJ, 954, 56
Olivares V., Picquenot A., Su Y., Gaspari M., Gendron-Marsolais M.-L.,
Polles F. L., Nulsen P., 2025, Nat. Astron., 9, 449

Osorno J., Nagar N., Richtler T., Humire P., Gebhardt K., Gultekin K., 2023,
A&A, 679, A37

Osterbrock D. E., Ferland G. J., 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei. University Science Books, Sausalito, CA

Owen F. N., Eilek J. A., Kassim N. E., 2000, ApJ, 543, 611

Peterson J. R., Fabian A. C., 2006, Phys. Rep., 427, 1

Pizzolato F., Soker N., 2005, ApJ, 632, 821

Polles F. L. et al., 2021, A&A, 651, A13

Proxauf B., Öttl S., Kimeswenger S., 2014, A&A, 561, A10

Reynolds C. S., Balbus S. A., Schekochihin A. A., 2015, ApJ, 815, 41

Rhea C. L. et al., 2025, AJ, 169, 203

Rich J. A., Kewley L. J., Dopita M. A., 2011, ApJ, 734, 87

Rich J. A., Kewley L. J., Dopita M. A., 2015, ApJS, 221, 28

Rousseau-Nepton L. et al., 2019, MNRAS, 489, 5530

Russell H. R. et al., 2017, MNRAS, 472, 4024

Salomé P., Combes F., Revaz Y., Edge A. C., Hatch N. A., Fabian A. C.,
Johnstone R. M., 2008, A&A, 484, 317

Sánchez-Blázquez P. et al., 2014, A&A, 570, A6

Sanders J. S. et al., 2016, MNRAS, 457, 82

Schure K. M., Kosenko D., Kaastra J. S., Keppens R., Vink J., 2009, A&A,
508, 751

Seon K.-I., Witt A. N., 2015, Highlights of Astronomy, 16, 576

Sharma P., McCourt M., Quataert E., Parrish I. J., 2012, MNRAS, 420,
3174

Simionescu A., Tremblay G., Werner N., Canning R. E. A., Allen S. W., Oonk
J. B. R., 2018, MNRAS, 475, 3004

Sparks W. B., Ford H. C., Kinney A. L., 1993, ApJ, 413, 531

Sparks W. B., Donahue M., Jordán A., Ferrarese L., Côté P., 2004, ApJ, 607,
294

Sparks W. B., Pringle J. E., Donahue M., Carswell R., Voit M., Cracraft M.,
Martin R. G., 2009, ApJ, 704, L20

Squire J., Chandran B. D. G., Meyrand R., 2020, ApJ, 891, L2

Sun M. et al., 2021, Nat. Astron., 6, 270

Tamhane P. et al., 2025, ApJ, 990, 36

Tremblay G. R. et al., 2018, ApJ, 865, 13

Vazdekis A., Sánchez-Blázquez P., Falcón-Barroso J., Cenarro A. J., Beasley
M. A., Cardiel N., Gorgas J., Peletier R. F., 2010, MNRAS, 404, 1639

Vigneron B. et al., 2024, ApJ, 962, 96

Voit G. M., 2018, ApJ, 868, 102

Voit G. M., Bryan G. L., O'Shea B. W., Donahue M., 2015, ApJ, 808, L30

Voit G. M., Meece G., Li Y., O'Shea B. W., Bryan G. L., Donahue M., 2017,
ApJ, 845, 80

Wang C., Ruszkowski M., Pfrommer C., Oh S. P., Yang H. Y. K., 2021,
MNRAS, 504, 898

Wang C., Oh S. P., Ruszkowski M., 2023, MNRAS, 519, 4408

Wei J.-J., Wu X.-F., Melia F., 2016, MNRAS, 463, 1144

Weil M. L., Bland-Hawthorn J., Malin D. F., 1997, ApJ, 490, 664

Werner N. et al., 2010, MNRAS, 407, 2063

Werner N. et al., 2013, ApJ, 767, 153

Yang H. Y. K., Reynolds C. S., 2016, ApJ, 829, 90

Young A. J., Wilson A. S., Mundell C. G., 2002, ApJ, 579, 560

Zhang C., Zhuravleva I., Gendron-Marsolais M.-L., Churazov E., Schekochihin A. A., Forman W. R., 2022, MNRAS, 517, 616

Zhuravleva I. et al., 2014, Nature, 515, 85

Zhuravleva I., Allen S. W., Mantz A., Werner N., 2018, ApJ, 865, 53

ZuHone J. A., Markevitch M., Zhuravleva I., 2016, ApJ, 817, 110

APPENDIX A: STELLAR CONTINUUM SPECTRA

Fig. A1 presents representative spectra and their corresponding stellar continuum fits used for continuum subtraction in the MEGARA and SITELLE data.
<figure><img src="image_14.png" /><figcaption>**Figure A1.** (a) Median spectrum of all spaxels in the MEGARA SE pointing, excluding edge regions affected by instrumental artefacts. Shaded areas indicate masked regions omitted from the PPXF fit (see the main text for details). (b) Integrated spectrum from a region within the central filament, extracted from SITELLE data at (α: 187.7075991; δ: 12.3938644)° with a circular aperture of radius ~1.4 arcsec. The spectrum shows data from the SN3 and SN1 filters, separated by a broken wavenumber axis but fitted simultaneously. In both panels, the best-fitting stellar continuum model from PPXF is shown as a thicker red line, with residuals plotted below. Main emission lines are marked by vertical black dotted lines, while major stellar absorption features are indicated in red dashed lines. In panel (a), the [O I]λ6300 line is visible, but too faint to be measured reliably on a spaxel-by-spaxel basis (SNR < 3); its line width appears consistent with the other emission features.</figcaption></figure>

## APPENDIX B: SITELLE VELOCITY DISPERSION MAP

Fig. B1 presents the velocity dispersion map derived from the SITELLE data (SN3 filter). We observe elevated values in the very

centre, reaching up to nearly 200 km s⁻¹, followed by intermediate dispersions of approximately 125 km s⁻¹ within the filaments covered by the MUSE WFM pointing. In the outer regions, including the FE filament, the dispersion drops to lower values of about ~80 km s⁻¹. However, we emphasize that these measurements are
<figure><img src="image_15.png" /></figure>

**Figure B1.** SITELLE velocity dispersion map showing spaxels with $SNR_{H\alpha} > 3$ and velocity uncertainties $\le 30$ km s$^{-1}$ ($2 \times 2$ binning). The MEGARA FOVs targeting the SE and FE filaments are outlined in blue boxes, and the FOV of the existing MUSE WFM observations is shown in a black dotted contour.

significantly affected by the limited spectral resolution of SITELLE, and should be interpreted with caution.

## APPENDIX C: IMPACT OF FILAMENT SEGMENTATION ON VSFS

To assess how the segmentation of filamentary substructures impacts the interpretation of turbulence, we performed a comparison of VSFs computed with and without applying the spatial segmentation described in Section 3.2.3. The goal is to test whether separating filaments leads to significant differences in the VSF parameters. Using the methodology described in Section 3.2.3, we computed localized VSFs across the MUSE velocity map with the same sampling (central positions spaced by ~100 pc and aperture radius of ~0.39 kpc), but without applying any subregion boundaries. For each central region, we then compared the resulting VSF to its counterpart obtained from the segmented approach by measuring the relative differences in three key parameters: the slope (power-law index), the flattening scale, and the amplitude at flattening.

Fig. C1 shows the variation in VSF parameter differences—slope, flattening scale, and amplitude—as a function of projected distance from the galaxy centre. The results indicate that the impact of segmentation on VSF measurements varies across the field. Differences in slope remain below 0.1 dex in approximately 80 per cent of the regions, while amplitude differences exceed $30$ km s$^{-1}$ in 28 per cent of the cases. Flattening scale differences greater than 0.1 kpc are observed in about 30 per cent of regions. Importantly, the most significant deviations in all three parameters tend to occur in the same areas—particularly in subregions F1, F4, and F5—where distinct kinematic components are known to overlap. This spatial coincidence suggests that the segmentation effectively reduces contamination from unrelated structures, leading to more reliable measurements of turbulence properties.

<figure><img src="image_16.png" /></figure>

**Figure C1.** (a) MUSE velocity map of the central filaments. Each square marks the central position of the localized VSF, computed within a circular aperture of radius 0.39 kpc, with and without applying the filament segmentation. White squares indicate locations where at least one VSF parameter exceeds the defined threshold (0.1 dex for the slope, $30$ km s$^{-1}$ for the amplitude, and 0.1 kpc for the flattening scale). (b) Variation of the VSF parameters as a function of projected distance from the galaxy centre. Horizontal lines mark the corresponding thresholds for each parameter. The vertical green line indicates the projected extent of the optical jet (~1.5 kpc).
## APPENDIX D: IMPACT OF SPAXEL SIZE ON VSFS

To assess whether differences in spaxel size between instruments impact the measurement of VSFs, we first compared SITELLE and MUSE data over matched subregions in central filaments of M87. SITELLE data have an original spaxel size of 0.32 arcsec × 0.32 arcsec (unbinned), while MUSE offers a finer spaxel size of 0.20 arcsec × 0.20 arcsec. The SITELLE unbinned map was produced following the same procedure described in Section 2.2, but without spatial binning. To enable a direct comparison at the same spaxel size, we also resampled the MUSE velocity maps to the SITELLE grid using the `reproject` function from the `ASTROPY` package. We manually defined five subregions (R1–R5) in the SITELLE velocity maps, selecting areas with well-defined, continuous filamentary structures, and sufficient pixel counts. The same spatial masks were then applied to the MUSE data. Owing to its finer spaxel size and higher SNR across the filaments, MUSE contains a larger number of spaxels per region. For each subregion, we computed the first-order VSF as described in Section 3.2. The results are shown in Fig. D1. Overall, the VSF slopes are broadly consistent between SITELLE and MUSE at intermediate and large scales. However, at small separations – particularly below the seeing limit of SITELLE (~1 arcsec) – SITELLE VSFs tend to diverge more strongly. Most regions show a consistent flattening scale between instruments, with some displaying a second flattening at larger separations where the number of pairs becomes insufficient. Across all subregions, SITELLE exhibits systematically higher VSF amplitudes. The resampled MUSE data produce nearly identical

VSFs to the original map, confirming that this amplitude offset is not due to spaxel size differences but likely stems from the noisier nature of the SITELLE velocity map, which shows enhanced pixel-to-pixel variations.

In a second test, we assessed the impact of spatial binning on the shape of the VSFs. For MEGARA, the binning was applied directly to the data cube by averaging the spectra within 2 × 2 and 3 × 3 spaxels, followed by refitting the binned spectra with a single-component model, as described in Section 2.1.2, to generate new velocity maps. For MUSE, where only the velocity map was available, binning was performed by averaging the velocity values within 3 × 3 and 10 × 10 pixel regions. VSFs were then recomputed for each binned velocity map using the same methodology as described previously. Fig. D2 presents the original and binned velocity maps, along with their corresponding VSFs. For the SE and FE filaments, the SITELLE VSFs (binning of 2 × 2) are also shown for reference. We find that binning has minimal effect on the overall shape of the VSFs. The flattening scale and its amplitude remain stable across all binning levels, including the extreme case of 10 × 10 binning in MUSE (effective spaxel size ~2 arcsec). The primary effect is a steepening of the slope at small scales, likely due to enhanced contrast between adjacent spaxels in the binned maps. Unlike Gaussian seeing, which smooths gradients, binning tends to preserve or amplify sharp transitions in the velocity field. These findings reinforce our conclusion that binning has moderate influence on VSF shape and amplitude. Once again, SITELLE systematically shows higher amplitudes than MUSE and MEGARA – even within comparable FOVs – further suggesting that the underlying data quality, rather than spatial sampling, drives this discrepancy.
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

**Figure D1.** Top panels: velocity maps of the central filaments obtained with SITELLE (left) and MUSE (right), limited to MUSE FOV. The subregions used for VSF analysis (R1–R5) are outlined in both maps. Bottom panels: VSFs computed for each subregion using SITELLE (solid lines) and MUSE (dashed lines). Reference slopes of 1/3 (Kolmogorov turbulence) and 2/3 are shown for comparison.
SE Filaments

a

Original
N = 778

<figure><img src="image_19.png" /></figure>

Binning 2 x 2
N = 221

<figure><img src="image_20.png" /></figure>

Binning 3 x 3
N = 114

<figure><img src="image_21.png" /></figure>

b

<figure><img src="image_22.png" /></figure>

FE Filament

a

Original
N = 531

<figure><img src="image_23.png" /></figure>

Binning 2 x 2
N = 152

<figure><img src="image_24.png" /></figure>

Binning 3 x 3
N = 76

<figure><img src="image_25.png" /></figure>

b

<figure><img src="image_26.png" /></figure>

Central Filaments

a

Original
N = 21 253

<figure><img src="image_27.png" /></figure>

Binning 3 x 3
N = 2694

<figure><img src="image_28.png" /></figure>

Binning 10 x 10
N = 330

<figure><img src="image_29.png" /></figure>

b

<figure><img src="image_30.png" /></figure>

**Figure D2.** (a) Velocity maps of the original and binned data sets for the SE filaments (top panel) and FE filaments (middle panel) with MEGARA, and central filaments (bottom panel) with MUSE. The number of spaxels *N* in each map is indicated in the upper left corner. Colour scale is consistent across maps. (b) Corresponding VSFs derived from each map, with reference slopes of 1/3 and 2/3 shown for comparison. Shaded areas indicate bins with an insufficient number of pairs in the original maps. For the SE and FE filaments, SITELLE VSFs (2 × 2) are also shown for comparison. The typical seeing limit (SITELLE: ~1.2 arcsec, and MUSE: 0.9 arcsec) is marked by a vertical dotted line.
## APPENDIX E: IMPACT OF THE WINDOW SIZE ON VSFS

To assess whether window size affects the measurement of VSFs, we performed a test using the full MUSE velocity map of the central region of M87. This approach is based on the subregions defined in Section 3.2.3. We defined a series of circular FOVs centred at representative locations, with diameter of 50, 75, and 100 spaxels (corresponding to 0.775, 1.1625, and 1.55 kpc) and computed the first-order VSF for each in their subregion using the method described in Section 3.2. The results are shown in Fig. E1.

We find that the overall shape of the VSFs does not show significant changes on smaller scales across different aperture sizes. Both the slope and the location of the first flattening are generally stable, indicating that small- and intermediate-scale turbulence is not strongly sensitive to our selected window size. However, as the aperture increases, the VSF extends to larger separations, and

a second rise and flattening becomes more pronounced – particularly in cases such as insets (a)/(b), where the first flattening gradually becomes less visible as large-scale motions begin to dominate. In other cases, such as inset (c), the first flattening becomes less distinct and is replaced by a smoother slope with increasing FOV.

This test confirms that, within the range of aperture sizes considered here, the measurement of VSFs is relatively robust up to scales of ~2 kpc. Nevertheless, larger FOVs increasingly may sample different dynamical regimes, shifting the dominant contribution from localized turbulence to more global kinematic structures. We also note that the window size can introduce artificial features in the VSF, typically at scales around half its length, although this effect is not always significantly apparent here. As such, care must be taken when interpreting the behaviour of the VSF or when comparing results across data sets with substantially different spatial coverage.
<figure><img src="image_31.png" /></figure>

**Figure E1.** The background shows the MUSE velocity map, with its scale bar located at the bottom right of the three circular apertures. Insets (a)–(e) display the VSFs computed within circular apertures of varying diameters – 50 (light-grey thick line), 75 (medium-grey solid line), and 100 spaxels (black dashed line) – centred on the black squares in a subregion (labelled F1–F5). These subregions are defined in the reference map at the top right. The relative size of each aperture with respect to the filament is illustrated by circles on the galaxy nucleus. Tick labels are displayed only in inset (d); all insets share the same axis scaling. Vertical blue lines mark half the size of each aperture’s spatial window, with line styles matching those of the corresponding circular apertures.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.