© The Authors 2026

1. Introduction

Galaxies with complex stellar kinematics have complex internal orbital structures. Predictions of stellar velocity maps (Statler 1991; Arnold et al. 1994) have indicated that triaxial galaxies could sustain diverse velocity structures because they can support multiple orbital families. In the separable potentials of the Stäckel type, there are four major orbital families: short-axis tubes, two types of long-axis tubes (inner and outer), and box orbits without net angular momentum (de Zeeuw 1985); these families have been visualised in Statler (1987)¹. In more realistic (non-rotating) triaxial potentials, irregular and resonant orbits, with descriptive names such as boxlets, bananas, saucers, fans, and fish, are also possible (Gerhard & Binney 1985; Schwarzschild 1993; Röttgers et al. 2014), and their spatial distributions are responsible for the observed shapes of galaxies (e.g. Schwarzschild 1979).

Non-axisymmetric galaxies can sustain velocity maps exhibiting kinematic twists, kinematically distinct cores (KDCs), and even counter-rotating cores. The origin of these kinematic structures was considered to be due to a combination of projection effects and the existence of different orbital families (Statler 1991; Arnold et al. 1994). However, except in a few cases where multiple slits were used to cover galaxies (e.g. Galletta 1987; Wagner 1990; Statler & Smecker-Hane 1999; Cappellari et al. 2002) and where velocity maps were obtained through interpolation between the slits, long-slit observations were not able to determine and characterise kinematic structures on the velocity maps.

The arrival of galaxy surveys with integral-field units (IFUs; Bacon & Monnet 2017) at the beginning of the 21st century enabled a more comprehensive view of the kinematics of early-type galaxies (ETGs; Emsellem et al. 2004). The IFU maps of the mean velocity, V, and the velocity dispersion, σ, covering up to one effective radius (R_e) have several important advantages over long-slit kinematics. On the theory side, the maps can be directly connected to the velocity dispersion tensor and the tensor virial theorem (Binney 2005; Cappellari et al. 2007). Observationally, they enable the characterisation of kinematic properties.

Krajnović et al. (2011) introduced five distinct classes based on a combination of visual inspection and harmonic analysis²: Group a, no detectable rotation (3%); Group b, irregular rotation (5%); Group c, KDCs (7%); Group d, counter-rotating discs (4%); and Group e, disc-like rotation (80%). The numbers in parenthesis refer to their frequencies in the magnitude-limited ATLAS^3D sample (Cappellari et al. 2011a) of local ETGs. Groups d and e are related, as their galaxies are made of stellar discs and are nearly axisymmetric, even when supporting counter-rotation. Galaxies in Groups a to c do not have stellar discs (Krajnović et al. 2013). Kinematic maps also allow for the calculation of the projected angular momentum as well as its proxy, the spin parameter ${\lambda }_{\mathrm{R}}\equiv ⟨R|V|⟩/R\sqrt{V^2+{\sigma }^2}$ (Emsellem et al. 2007), which can be calibrated to distinguish between galaxies with and without stellar discs (Krajnović et al. 2011; Emsellem et al. 2011; Cappellari 2016; van de Sande et al. 2021b).

Numerical simulations (Röttgers et al. 2014; Frigo et al. 2019) and dynamical models constrained by IFU data (van den Bosch et al. 2008; den Brok et al. 2021) are able to connect the observed complexity of stellar kinematics with the distribution of internal orbits. They show that velocity maps should be understood as luminosity-weighted superpositions of stellar motions stemming from different orbital families, including those with no net streaming. Orbital families can be co-spatial, with stars moving in the prograde and retrograde directions within any family, cancelling or reinforcing streaming observed in the velocity maps. In some cases, velocity maps can appear as if there are multiple spin reversals (den Brok et al. 2021). Crucially, the complexity of the observed kinematics can be related to the merger history (Naab et al. 2014; Ebrová & Łokas 2017; Li et al. 2018; Ebrová et al. 2021; Lagos et al. 2022).

In this work, we present the MUSE Most Massive Galaxies (M3G) Survey consisting of galaxies brighter than −25.5 in the K band. We focus on the survey design and observations, presenting the kinematics and an analysis of the velocity maps. Such galaxies are rare in the nearby Universe, and their stellar kinematics are not yet well mapped, in spite of recent works (Loubser et al. 2008; Ma et al. 2014; Loubser et al. 2022). An expectation at the time of the M3G survey definition was that they would have little or no rotation, thus belonging to Groups a - c. If there were any exceptional kinematic structures, they would be confined to their cores, within a few kiloparsecs at most. The primary purpose of this paper is to demonstrate and quantify the surprising richness of observed features in M3G velocity maps. As a result, the richness of the kinematic features confirms their origin through the superposition of stellar motions belonging to different orbital families.

This paper is organised as follows: In Section 2 we present the M3G sample, MUSE observations, and the data reduction. In Section 3 we describe the extraction of stellar kinematics, while in Section 4 we present the kinematic analysis and the main results of the paper. A discussion and the conclusions are presented in Sections 5 and 6, respectively. The kinematic maps of M3G galaxies and other supporting material are attached in the Appendices. We used the standard cosmology and H₀ = 70 km/s/Mpc.

2. The galaxy sample, MUSE observations, and data reduction

The M3G sample, MUSE observations and the data reduction were already partially presented in Krajnović et al. (2018). Here we repeat the pertinent details and point to the new developments since the publication of that paper. Some of these developments were already discussed in Pagotto et al. (2021) and den Brok et al. (2021, 2024) and are summarised here.

2.1. The properties of M3G sample galaxies

The idea behind the M3G sample was to explore the spectroscopic properties of galaxies in a way not possible by previous surveys: to observe very massive galaxies that are rare in the nearby Universe beyond their effective radius. The sample was in practice selected by imposing the following criteria:

• Galaxy sizes were required to match the MUSE field of view (FoV) such that 2R_e are covered within a single pointing.

• Galaxies were required to be brighter than −25.5 magnitude, corresponding to ∼5 × 10¹¹ M_⊙ (following Eq. (2) in Cappellari 2013), the upper limit brightness and mass ranges for the ATLAS^3D sample (Cappellari et al. 2011a).

• Galaxies were supposed to be in dense environments and, in particular, rich clusters (richer than the Virgo cluster, also covered by the ATLAS^3D sample).

• Most galaxies of this brightness limit were expected to be the brightest cluster galaxies (BCGs), but a significant fraction of galaxies was required to belong to non-BCGs.

• The spatial resolution was required to be better than 1 kpc, as this is the typical size of KDCs in massive ETGs (Krajnović et al. 2011). Assuming natural seeing conditions at the VLT, this constraint limited the distance to about 250 Mpc for the selection of massive galaxies.

We searched for candidate galaxies in the core of the Shapley Supercluster (SSC; Shapley 1930; Merluzzi et al. 2010, 2015; Haines et al. 2018), focusing on Abell 3556, 3558 and 3562, the three largest clusters within the overall supercluster structure. The magnitude selection provided 14 galaxies, among which are also BCGs of the three clusters. We further matched this SSC sub-sample with 11 BCGs based on the sample of galaxies presented in Laine et al. (2003), which itself is based on the survey of Abell et al. (1989) clusters by Lauer & Postman (1994) and Postman & Lauer (1995). The galaxies were located in the southern hemisphere (declination < 5°), in clusters that have the Abell (1958) richness³ larger than the Virgo Cluster (33), that match the magnitude and distance cut of the SSC subsample. As emphasised by Abell et al. (1989), richness is not a substitute for local luminosity functions, and we use it only as a relative measure of environments denser than those covered by the ATLAS^3D survey (Cappellari et al. 2011b). The three BCGs in the SSC sample are also part of the Laine et al. (2003) sample, and therefore all 14 BCGs have HST imaging, while only two non-BCGs were imaged with HST.

The final M3G sample is made of 25 galaxies, 14 BCGs and 11 members of the three main galaxy clusters in the core of the SSC. We refer to these galaxies as non-BCGs, and they comprise at least the second and the third brightest galaxies in their respective clusters, as well as galaxies ranked to lower brightnesses. Based on galaxies in nearby clusters (e.g. Virgo or Coma cluster), the second brightest galaxies are often very similar to the nominal BCGs, and for some interpretations this needs to be kept in mind (Lauer et al. 2014). The general properties of the sample galaxies are presented in Table 1. The emission-line properties and the dust content were discussed in Pagotto et al. (2021), while the stellar population and elemental abundances were presented in den Brok et al. (2024). We do not repeat them here.

Three galaxies in the M3G sample have incorrect designation in previous papers (Krajnović et al. 2018; Pagotto et al. 2021; den Brok et al. 2024). They are PGC 097958, PGC 099188 and PGC 099522. The correct designations for PGC 097958 and PGC 099522 are LEDA 097958 (in Abell 3558) and LEDA 099522 (in Abell 3562), because the Catalogue of Principle Galaxies has 73197 official entries (Paturel et al. 1989), while catalogue entries with larger number are designated as LEDA. They are listed in the Lyon-Medon Extragalactic Database (Paturel et al. 1988), and services such as SIMBAD⁴ easily resolve between PGC and LEDA catalogue names. The galaxy previously named PGC 099188 is actually PGC 047548 (in Abell 3562), with coordinates RA = 13 31 27.52 Dec = −31 49 14.60 (J2000), also designated as ShaSS403013239 (Mercurio et al. 2015). PGC 099188 (LEDA 099188), or ShaSS403010664 is a faint galaxy at z = 0.21 (Haines et al. 2018). In order to highlight this inaccuracy, but to simplify the cross-correlation with previous M3G sample papers and designation of the MUSE dataset in the ESO archive⁵ we nevertheless keep the designation of PGC099188* (highlighted by a star) for PGC 047548.

We made use of the 2MASS K_s band magnitudes (k_m_ext, the extrapolated K_s band 2MASS magnitude) and size measurements from the All Sky Extended Source Catalogue (XSC; Jarrett et al. 2000; Skrutskie et al. 2006). The R_e of our galaxies was estimated as the mean of the 2MASS measured effective radii in J, H and K_s bands, R_e = mean(j_r_eff, h_r_eff, k_r_eff). We compared these R_e to R_e^maj obtained as the major axis of the half-light isophote (Cappellari et al. 2013), by using the uniform r band VST images of galaxies in the SSC (Merluzzi et al. 2015) to construct multi-gaussian expansion (MGE) models⁶ (Emsellem et al. 1994), and estimate R_e^maj. Comparing the calculated R_e^maj with the R_e for the 14 SSC galaxies, we obtained a linear relation, with the slope b = 1.0 ± 0.1. As there is no uniform VST imaging for the full sample to facilitate a uniform calculation of R_e^maj, we adopted the 2MASS-based half-light radii, as defined above.

A major source of uncertainty for the sample properties is the distance measurements. Only about 1/3 of the sample has known distance estimates, and those are based on scaling relations, such as the fundamental plane. Therefore, we estimated the distance as follows. Firstly, we collected all available distances from the NASA Extragalactic database⁷. When there were two or more distance estimates, which were not strongly discrepant, we used the average value. For galaxies that did not have distance estimates, we collected the redshift measurements and used the Hubble law to obtain the distance given the redshift:

$$\begin{array}{c}\hfill D\simeq \frac{c}{H_0}\frac{c{(z+1)}^2-1}{c{(z+1)}^2+1}.\end{array}$$(1)

Finally, we required that all galaxies within the same cluster (SSC galaxies) have the same distance, ignoring the distribution of galaxies within a cluster. We calculated the SSC clusters distances as the median of the available distances of all cluster galaxies present in our sample. Based on the mean recession velocities (Haines et al. 2018), the distances to Abell 3556, Abell 3558 and Abell 3562 are 200 Mpc, 202 Mpc and 205 Mpc, which is up to 10% different from our final estimates (Table 1). Haines et al. (2018) report a large spread of velocities for the three SSC clusters (σ_c= 628 km/s, 1007 km/s and 769 km/s, respectively), and a distance range of about 20–30 Mpc within the cluster, also amounting to about 10% with respect to the cluster mean distance. Overall, we estimate the accuracy of the distances to be about 10%, but this has no influence on the results of this paper. The distribution of distances is shown in Fig. 1: SSC galaxies are located between 180 and 200 Mpc, while the BCGs have a somewhat larger spread. The adopted distances mean that 1″ is between 0.73 and 1.12 kpc for the M3G sample. The new distances and the effective radii take precedence over those published in previous M3G publications.

Fig. 1. Properties of the M3G sample. Top left: Histogram of the distances to the BCGs (orange) and the satellites galaxies (blue). Bottom left: Richness of the host galaxy cluster versus the absolute Ks magnitude of the corresponding BCG. BCGs in the SSC are shown as pentagrams (Abell 3558, Abell 3556 and Abell 3562 in order of decreasing brightness of the BCGs). The Virgo and Coma cluster (squares) have been added as references for the ATLAS3D and the MASSIVE Surveys. Right: Size–luminosity relation for the M3G sample galaxies, split into BCGs (red) and satellites (blue). Galaxies from two other surveys have been added for comparison: ETGs from the ATLAS3D survey (light green large squares), spirals from the ATLAS3D survey parent sample (dark green small squares), and ETGs from the MASSIVE survey (brown triangles). ATLAS3D sample is a magnitude-limited sample within 42 Mpc, while the galaxies in the MASSIVE are the brightest galaxies within 100 Mpc. The solid and dashed black lines are relations for fast rotators and S0-Sa galaxies from the ATLAS3D survey, respectively. The red line denotes the so-called ‘zone of avoidance’. All three relations are taken from Cappellari et al. (2011a).

Figure 1 shows where the M3G targets stand with respect to local samples of ETGs investigated with IFU spectroscopy. M3G galaxies are found on the apex of the luminosity-size relation, reaching luminosities of 10¹² L_⊙ and half-light sizes larger than 20 kpc. They are thus similar to the galaxies in the MASSIVE sample (Ma et al. 2014), with a few BCGs being somewhat more luminous. The non-BCGs in the sample have, as expected, lower luminosities and smaller sizes, but can still be considered extreme systems.

In terms of the environment, M3G galaxies are found in clusters with richnesses between those of the Virgo and Coma clusters. Abell 3558 is by far the richest cluster, containing almost a factor of two more galaxies than other probed clusters, and contains the brightest galaxy in the sample (PGC 047202). Abell 3562, while being the second richest cluster in the sample, harbours the faintest BCGs (PGC 047752). Two other galaxies from this cluster, PGC 099188* and PGC 099522, were found to actually have larger 2MASS apparent luminosities than PGC 047752 (10.7 mag and 10.6 mag, compared to 10.8 mag, respectively). This could be driven by the presence of nearby overlapping satellites: they have significantly different systemic velocities, but fall within the MUSE FoV and strongly influence the light distributions within the central 10″. Additionally, Abell 3562 experienced a recent supersonic encounter with the SC 1329-313 group (Finoguenov et al. 2004). This is relevant as the centre of the X-ray emission in the vicinity of SC 1329-313 group coincides with the location of PGC 47548, while the similarly bright PGC 047590 is located in the vicinity along the same X-ray filament emission. A brighter X-ray emission source is centred on PGC 047752 identified as the BCG of Abell 3562 by Postman & Lauer (1995). Haines et al. (2018) assigns both PGC 47548 and PGC 047590 to Abell 3562 based on velocity-related arguments. Therefore, M3G members of Abell 3562 have similar apparent 2MASS K_s band magnitudes (10.7, 10.6, 10.8 and 10.6 mag in the order of the increasing PGC number), where the faintest can be taken as the local BCG, while the brightness of others might be overestimated by the presence of companions.

We note that the differences in the derived absolute magnitudes for M3G Abell 3562 galaxies are less than 0.3 mag, a level sometimes used to consider the second-ranked galaxy a close rival to the local BCG (Lauer et al. 2014). In this respect, all M3G galaxies of Abell 3562 could be seen as close rivals for the BCG place. We nevertheless keep PGC 047752 (also known as ESO 444-G72) as the BCG of Abell 3562. PGC 099188*, PGC 047752 and PGC 099522 are assigned the rank of the second, third and fourth galaxies in the cluster, and we remind the reader to keep in mind the connection with SC 1329-313 group for PGC 099188* and PGC 047590, and a possibly overestimated luminosity for PGC 099522.

2.2. MUSE observations

The M3G sample was observed as part of the MUSE Guaranteed Time Observations (GTO) during ESO periods 94-99 spanning three years (2014–2017)⁸. Each galaxy was observed using MUSE in Wide Field Mode without Adaptive Optics and a nominal wavelength range (480–930 nm). Observing Blocks (OB) were split into four on-target and two sky observations within an OSOOSO sequence (O is Object, S is Sky). Consecutive on-target observations were rotated by 90^o and dithered to reduce systematics between the individual spectrographs. The sky exposures (120s) were pointing several arc minutes away from the galaxy to areas free from cluster galaxies (and foreground stars). No other special calibration exposures were used, except routine standard stars observed once per night.

The M3G survey was assigned to be a “poor seeing” programme within the GTO projects, meaning that it was executed when the FWHM of the ambient seeing would exceed 1″, but be no worse than $1\stackrel{″}{.}4$. A fraction of GTO time was further allocated for relatively short on-target exposures (≈1800 s) using “good seeing” (GS: FWHM < $0\stackrel{″}{.}8$) conditions to resolve the central core structures (smaller than 1 kpc). In practice, exposures within one OB exhibit strongly varying seeing conditions. In Fig. 2, we show the distribution of recorded seeing values by the VLT auto-guider. The resulting mean seeing of the survey is 0.93″, with 28% of the sample galaxies having the mean seeing poorer than 1″. This is better than expected given the observational strategy. In several cases, however, the GS exposures did not achieve the target specification, with 44% galaxies having seeing worse than the targeted $0\stackrel{″}{.}8$, but the median seeing of the GS data cubes is $0\stackrel{″}{.}75$. The seeing values for each galaxy are listed in Table 2. In this work, we combine all available exposures to increase the signal and only occasionally use the GS cubes to verify our kinematic analysis in the nuclei (Section 4.1.2), but do not present them here.

Fig. 2. Distribution of the VLT auto-guider seeing measurements averaged over each exposure taken with MUSE for all M3G galaxies. Horizontal lines denote the two seeing limits imposed on the programme. Blue squares show the mean value of all observations for each galaxy. Red circles are exposures targeted to be done in ‘good seeing’ conditions (FWHM < 0.8″).

Total exposure times for each galaxy ranged from 2 to 6 h. Table 2 also lists the limiting surface brightness magnitudes reached at the edge of the extracted kinematic maps (see Section 3). The magnitudes were extracted from the MUSE reconstructed images using the SDSS r band filter. The average limiting surface brightness reached is 23.21 mag/″² (AB magnitudes) with a standard deviation of 0.38 mag/″², showcasing the achieved uniformity of observations.

Each galaxy was observed by placing its nucleus at the centre of the MUSE FoV. For the central galaxy of Abell 3558 and the brightest in the sample, PGC 047202, we built a 2 × 2 mosaic from the nominal “poor-seeing” exposures, each time placing its nucleus in one MUSE corner, in addition to “good-seeing” exposures with the nucleus placed in the centre on the MUSE FoV.

The quality of the individual exposures is typically high. A few individual exposures were not completed or, upon inspection, found to be of insufficient quality to be included for the final scientific analysis. In some cases, exposures had a reduced flux (i.e. down by 20–40% for PGC 047202 in the SW quadrant, see Fig. C.1), due to passing clouds, and in two instances for PGC 065588, satellite trails were also recorded. These observations were nevertheless kept and combined with the rest, resulting in still valuable albeit not ideal data. Exposures were rejected from the final creation of the data cubes if they had been aborted, usually due to the loss of a guide star (PGC 03342, PGC 15524, PGC 047202). One OB belonging to PGC 048896 was observed through clouds and resulted in no exploitable data, while two PGC 019085 OBs suffered from incorrect offsets when moving between on-target and sky locations.

2.3. Data reduction

The M3G galaxies were observed during 15 observing runs spread over six semesters, with data reduction proceeding on individual exposures as they were completed. This means that we used different versions (v1.2 to v1.7) of the MUSE data reduction pipeline Weilbacher et al. (2020) as they became available. We verified that this has not significantly impacted the reduction products and the results presented in this paper. All reductions followed the standard MUSE pipeline steps, from master calibration files (bias, flat fields, trace tables), including the wavelength calibration and the line-spread function (LSF) estimates for each spectral layer. In a number of cases, we used twilight flats (if observed on the same night). Instrument geometry and astrometry were produced at each observing run by the GTO team. We paid particular attention to using the illumination flats as close to the individual on-target exposures as possible.

Removing the sky contribution is non-trivial when objects fill the full MUSE FoV. We used the dedicated sky observations to build the spectra consisting of sky lines and sky continuum. We applied the sky subtraction method “model” within the scientific post-processing pipeline stage, which re-fits the sky lines before subtracting them from the on-target data. The sky continuum is subtracted directly without further modifications. Sky spectra were always associated with the closest in time on-target exposures. Such sky subtraction worked well in the blue part of the spectra (< 700 nm) and is sufficient for this study. Nevertheless, the process was less successful in the red part, and we point the reader to den Brok et al. (2024) for further details.

The flux calibration was calculated using the ESO-approved MUSE standard stars, which were in almost all cases taken at the same night as the M3G observations, and reduced with the MUSE pipeline recipes. As we did not observe telluric stars, our initial (Krajnović et al. 2018) telluric correction was based on the standard stars and was not effective. We improved the telluric correction by using the MOLECFIT software (v1.5.9; Smette et al. 2015), which computes a theoretical absorption model based on radiative transfer and atmospheric molecular line database. We first ran the scientific post-processing part of the pipeline with telluric correction switched off, preparing data cubes with spectra containing the telluric absorption features. We selected spectra of foreground stars, as featureless as possible, as inputs for MOLECFIT. The spectra were taken from merged data cubes combining all on-target exposures (see below). A few galaxies did not have suitable stars, and we used instead the central spectrum of the galaxy as the input for MOLECFIT to calculate the telluric absorption. These telluric spectra were then applied as a correction to all spectra in the merged data cubes. As shown in Pagotto et al. (2021) and den Brok et al. (2024), this procedure significantly improves the spectra and enables the analysis of emission lines and chemical properties of M3G galaxies. In this work, we still track the location of the strong sky and telluric features and mask them during the extraction of stellar kinematics.

The number of exposures per galaxy that were deemed of sufficient quality to be combined is listed in Table 2. We created “white” light images from each data cube by summing the cube along the wavelength dimension, and used them to provide the relative spatial offsets between the exposures. These offsets (and relative flux weights) were passed to the pipeline recipe for merging (MUSE_EXP_COMBINE), which was also used to reject the cosmic rays, applying the ‘median’ rejection method. For PGC 047202, where we combined 101 images, over the area of almost 2″ × 2″, the merging required a somewhat different approach. We determined the offsets between the individual exposures as before and defined a World Coordinate System with the location of the four quadrants of the mosaic. Then we used MUSE_EXP_COMBINE to combine all 101 exposures, but split into seven wavelength regions, each 65 nm long. This resulted in seven combined data cubes with short wavelength ranges but full spatial extent. Those were further combined into a final cube covering the full spectral range, using MUSE_CUBE_COMBINE⁹.

We note that until this point, all calibrations and corrections (including the sky subtraction) were done on the so-called MUSE Pixel Tables. The product of the merging routines is a combined data cube, which we use from this point for further analysis. The telluric correction was, however, not done on the pixel tables, but on the merged (non-telluric-subtracted) cubes. The final data cubes of the M3G sample are oriented with North up and East to the left, have 0.2″ × 0.2″ spatial sampling, and a spectral sampling of 1.25 Å.

3. Stellar kinematics extraction

The M3G galaxies are in dense environments and often surrounded by smaller satellites, which partially obstruct the view of the main galaxy. We have mitigated the impact by performing several different extractions of the stellar kinematics via the penalised Pixel Fitting method¹⁰ (pPXF; Cappellari & Emsellem 2004; Cappellari 2017). These included a pPXF fit to a global spectrum, extracted within an elliptical “effective” aperture defined by the R_e and the galaxy flattening (Table 1). This was followed first by a single component fit, and then by a multiple component fit to each spectrum, with two different parametrisation of the line-of-sight velocity distribution (LOSVD). The extraction presented here is broadly similar to that of Krajnović et al. (2018), but it differs in several key points, superseding the previously published kinematics. Nevertheless, the main result of Krajnović et al. (2018), pertaining to the distribution of the kinematic misalignment angles, remains valid.

3.1. Spatial binning and pPXF setup

Starting from the final telluric-corrected data cubes, we proceeded to spatially bin them to obtain a uniform signal-to-noise ratio (S/N) across the FoV. We first estimated the S/N for each spectrum in the cube using regions of the spectra lacking emission lines, sky or telluric residuals and strong absorption lines, with the noise estimated from the pipeline-produced variance spectra. Contrary to the previous approach, we only spatially masked prominent stars or compact sources, but included the satellite galaxies. Stars and compact objects to be masked were automatically detected via the subtraction of a smooth MGE model of the “white light” images: associated masks were built using a 1 − σ threshold above the median residual value. We also excluded all spectra that had an estimated S/N < 2: those are located at the edges of the MUSE FoV and are dominated by sky subtraction residuals. We used the Voronoi binning method¹¹ (Cappellari & Copin 2003), and set the target S/N to 50.

For galaxies that have no satellites, the pPXF setup was relatively standard: we parametrised the LOSVD with Gauss-Hermite moments (Gerhard 1993; van der Marel & Franx 1993), with the mean velocity V, the velocity dispersion σ describing the Gaussian part, and the Hermite coefficients h₃ and h₄, describing the asymmetric and symmetric departures from the Gaussian, respectively. We used additive polynomials of the 12th order, masked all emission lines and a few spectral regions likely to have sky or telluric residuals. Before fitting, the spectra were de-redshifted using an estimated redshift (equivalent to typically ≈20 nm for galaxies in our sample), and the fit was constrained between 470–680 nm (rest frame). We used the MILES stellar library (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011) as stellar templates, convolved to follow the MUSE line-spread function (LSF; Guérou et al. 2017). As a first step, we fitted the sum of all spectra within the elliptical effective aperture with the full set of MILES stellar templates, producing an optimal template. This optimal template was used in the subsequent pPXF run to fit all individual Voronoi bins, considerably speeding the fitting process and error estimation. We performed a test by generation optimal templates for each bin, but the resulting kinematics was not different compared to the nominal. Furthermore, as a test we extracted kinematics of a few galaxies using the higher spectral resolution X-shooter Stellar Library (Verro et al. 2022), which produced essentially identical kinematics and we opted to keep using the MILES templates. Finally, we tested extraction using the Calcium Triplet complex (typically redshifted to about 900 nm), but the extracted kinematics did not improve compared to the nominal, and the velocity moments had larger errors due to the increase noise caused by the residuals of the telluric correction.

The presence of overlapping satellites made the extraction of the stellar kinematics quite challenging, as some bins could present a mixed contribution of both the main target and a satellite. We took advantage of the different systemic velocities of those two contributions to implement a multi-component fitting process. Our approach is similar to the one commonly used for separating kinematic components within a given target (e.g. Coccato et al. 2011; Johnston et al. 2013; Coccato et al. 2018), except for the fact that we do not use any priors on stellar populations or light distributions (e.g. Johnston et al. 2012; Schmidt et al. 2019). When fitting multiple components, we kept the same pPXF set-up as outlined above, and allowed for one ‘main’ component with a Gauss-Hermite LOSVD and other components (satellite) described by a single Gaussian LOSVD. Optimal templates were built separately for the central object and the satellites using high S/N regions for each of those. These fits further provide the starting kinematic values (V and σ) that were used as a constraint for the subsequent multiple-component runs.

Once both the single and multiple components fits were done for each spatial bin of the MUSE FoV, we looked at the reduced χ² statistics of the pPXF results to determine which spatial bin belongs to the main galaxy (i.e. dominated by its light), to an intruding satellite, or its spectra could be decomposed into two components. We did this by imposing σ_sat ≥ σ_inst ≈ 50 km/s in each bin and the following restrictions on χ_R² = χ_single²/χ_multi², the ratio of the χ² obtained from the single and multi-component pPXF fits. Bins with χ_R² < 1.03 we assigned solely to the main galaxy. Those with 1.03 < χ_R² < 1.05 were deemed to carry both contributions from the galaxy and a satellite which could be disentangled. Our tests showed that bins with χ_R² > 1.05 provide unreliable kinematics for the main galaxy. These bins are dominated by the satellite’s light, their kinematics are typically robust, and we assigned them to the satellite.

In Fig. 3 we show an example pPXF fit to three different spectra coming from PGC 099188*, which has a large interloper satellite galaxy in the north east region. The top and bottom panel spectra belong to the main galaxy and the satellite, respectively, while the middle panel spectrum is taken from a location approximately midway between the centres of the galaxy and the satellite. The spectral features in the middle panel clearly indicate that the spectrum is composed of two components, which can be separated into a galaxy spectrum and a satellite spectrum. The interloping satellite is approximately 2000 km/s ahead of PGC 099188*, which is quite similar to some other cases in the sample. The main galaxy features are, however, not visible in the spectrum taken from the central region of the satellite. Nevertheless, the multi-component pPXF fit can recover the kinematics of the sample galaxy, which is of interest for our study, over a large spatial range. The middle panels of Fig. S2 in the supplementary material show the results of the multi-component fit in the case of PGC 099188, and the separation of the kinematics of the main galaxy and the satellite. More details about the separation of the satellite kinematics are provided in Appendix A.

Fig. 3. Spectral fits using pPXF for three different spectra from the PGC 099188* MUSE field, which features a large interloping satellite. From top to bottom, the panels show the central spectrum of PGC 099188*, a spectrum from the midpoint between the centre of PGC 099188* and the interloping satellite, and the central spectrum of the interloping satellite. In each panel, the black line shows the observed spectrum, and the red line shows the pPXF fit, with the residual of the fit in green points at the bottom. Green shaded polygons are masked regions excluded from the fit. The middle panel has two additional spectra shown by light and dark blue colours, belonging to the main galaxy and the satellite, respectively. These are model spectra vertically shifted for an arbitrary amount, obtained from the pPXF fits highlighting the satellite and the galaxy contributions to the total spectrum.

Preliminary mean velocity maps of the M3G galaxies were presented in Krajnović et al. (2018). Here we present the improved mean velocity maps, together with maps of velocity dispersions and h₃ and h₄ Gauss-Hermite moments. All 25 galaxies can be seen in Fig. C.1. In a few cases, there are notable differences between the new and old velocity maps, mostly in the regions covered by the outer bins and typically beyond 1R_e: those are mostly driven by an improved accounting for the zeroth order velocity term (Section 4.1.1). We note the measurements of the global velocity dispersions σ_e, presented in Table 1, were done by summing all spectra within a half-light aperture and parametrising the LOSVD with a Gaussian distribution, applying the single-component pPXF set-up.

4. Kinematic analysis

In this section, we analyse different aspects of the extracted M3G kinematics. We start with the ‘kinemetric’ analysis to describe and classify the features on velocity maps. We measure the spin parameter as a proxy for stellar angular momentum and classify the M3G galaxies into slow and fast rotators. Finally, we look into the properties of the third Gauss-Hermite moment and investigate the origin of regular rotation in massive galaxies.

4.1. Kinemetric analysis of velocity maps

We applied the KINEMETRY method (Krajnović et al. 2006) stepping through a velocity map at different radii and defining an ellipse at each radius along which the azimuthal variation of velocities are best reproduced by a simple cosine V(R) = V_rot(R)×cos(θ). The optimisation was done via harmonic analysis following

$$\begin{array}{c}\hfill V(R,\theta )=A_0(R)+{\displaystyle \sum _{n=1}^N}{A}_n(R)sin(n\theta )+{\displaystyle \sum _{n=1}^N}{B}_n(R)cos(n\theta ).\end{array}$$(2)

For odd moments, such as the mean velocity map, B₁(R) = V_rot(R) = V_csin(i). Therefore, harmonic terms up to N = 3, with the exception of B₁, were fitted to define the parameters of the ellipse, while the N ≥ 5 terms were not fitted but instead calculated once the best-fitting ellipse had been defined. The higher order terms represent the deviation from the assumed cosine law azimuthal velocity variations. Further, A₀ is a constant term, and it can be subtracted from each ellipse ring to reveal the relative motions within the galaxy.

When analysing maps of stellar velocities, it is useful to define amplitude terms as $k_n=\sqrt{A_n^2+B_n^2}$ (Krajnović et al. 2008). The term k₁ is dominated by B₁ and describes the rotation, while k₅, which is based on the first pair of harmonic terms that are not minimised, is typically used to estimate the deviation from the cosine law. Krajnović et al. (2011) required that a unit-less quantity $\overline{k_5/k_1}\le 0.04$ for a velocity map be well described by the cosine law, where the average is within the effective radius. Such maps are characterised as having ‘regular rotation’, while maps with $\overline{k_5/k_1}>0.04$ are described as having ‘non-regular rotation’. We also note that in this notation, k₀ = A₀.

We ran KINEMETRY in several setups and on two types of maps tracing the LOSVD moments: on the flux images (performing classical isophotometric analysis; Lauer 1985; Jedrzejewski 1987) and on the mean velocity maps. The analysis of the velocity maps were separated into two steps, the first one dealing with the removal of the zeroth kinemetry term (related to the systemic velocity) and the second with the actual harmonic analysis of the corrected velocity maps.

4.1.1. Spatial variations of the k₀ term

The outer regions of some velocity maps extracted in Section 3 display an unusual symmetric distribution of velocities with respect to the principal axes. Two examples are presented in the first column of Fig. B.1. Velocity maps of massive galaxies are expected to exhibit no rotation across the field (Group a), irregular and low amplitude rotation (Group b), or spatially localised rotating components exhibiting point-antisymmetric velocities fields (Group c). The velocity maps in Fig. B.1, however, have both non-zero velocities and are ‘symmetric’ at large radii. KINEMETRY analysis indicates that in these galaxies the k₀ term is changing with radius, independently from the higher terms, notably k₁, which traces the rotation around the galaxy centre.

A similar symmetric change in velocities was also noticed by Bender et al. (2015) in NGC 6166, the BCG of Abell 2199, where the velocities at larger radii approach the average velocity of the cluster galaxies. A possible explanation is that the halo of the BCG consists of tidal debris, which is more at rest with respect to the cluster than the central galaxy. If we combine the systemic velocities¹² of galaxies on Fig. B.1 with their k₀ values, V_sys + k₀ approach the values of the respective mean cluster velocities. For example, PGC047202 V_sys is 700 km/s below the cluster mean velocity V_c = 14450 km/s (Table S2 of the supplementary material), still within the spread of the cluster velocities σ_c = 1007 km/s (Lauer et al. 2014). At the edge of the MUSE field, even though only at 2R_e and well within the main galaxy, there is a significant decrease in the difference to ∼590 km/s.

The majority of M3G galaxies has nearly constant k₀ across the field. We quantify the absolute radial change as δk₀ = |k₀(0)−k₀(R)|, and introduce the maximum change Δk₀ = max(k₀(R)) − min(k₀(R)) as the difference between the maximum and the minimum k₀ values¹³. We show the radial change δk₀ in Fig. 4 and list Δk₀ in Table 1. There are seven galaxies that have Δk₀ > 20 km/s (28%), five of which are BCGs and two are non-BCGs.

Fig. 4. Radial profiles of the absolute k0 variation for M3G galaxies (δk0 = |k0(0)−k0(R)|) compared to the same from the ATLAS3D galaxies. The M3G galaxies are divided into BCGs (red lines) and non-BCGs (blue lines), while ATLAS3D galaxies are divided into fast (light blue dashed lines) and slow (light red dashed lines) rotators. The regions shaded purple show the standard deviation of all ATLAS3D and M3G (darker shaded regions) velocity errors. Names highlight specific galaxies with the largest difference between the maximum and the minimum values of k0: Δk0 > 20 km/s.

We compare in Figure 4δk₀ with those from the ATLAS^3D sample from Krajnović et al. (2011). The M3G galaxies are separated into BCGs and non-BCGs, while the ATLAS^3D are separated into fast and slow rotators. The majority of galaxies in both samples do not show large variations in δk₀. Only 19 ATLAS^3D are found to have Δk₀ > 20 km/s (7% of the sample).

In the ATLAS^3D sample notable outliers are NGC 4486 and NGC 4753 (in order of decreasing Δk₀). NGC 4486 is a massive slow rotator with an active SMBH and complex stellar kinematics (Emsellem et al. 2014). δk₀ is relatively constant at ∼35 km/s, except for a sudden jump between the central value and the values at radii larger than 2″. The nucleus of NGC4486 is dominated by AGN activity (for examples of central spectra of NGC4486 at different spatial resolutions and the influence of the AGN emission see Simon et al. 2024). Given that this type of δk₀ profiles are rare, the issue is likely related to difficulties in extracting reliable stellar kinematics from the nucleus. On the other hand, NGC 4753 is a galaxy with disturbed large scale morphology, indicating a recent merger (Bílek et al. 2020), and its δk₀ changes gradually with radius. NGC 4486 and NGC 4753 could be taken as showcases for two different reasons of k₀ radial variations: due to difficulties in extracting the correct nuclear velocities¹⁴ and due to disturbed equilibrium at large radii.

Figure 4 highlights seven M3G galaxies of which the largest δk₀ variation is seen for PGC 047202, the brightest and the largest galaxy in the sample, also found in the densest environment. A counter-example is PGC 073000, characterised with a sudden jump within the central 2″ reaching a value similar to NGC4486, which stays constant within ∼10″ and then steadily declines. The central sudden variation could be related to the contamination of the nuclear spectra, as this galaxy has very complex gas kinematics, ionised by an AGN (Pagotto et al. 2021). PGC 015524, with the second highest change in V_sys, has a similar behaviour to PGC 047202, with a steady increase of δk₀ from the central value. The increase is negligible (less than the typical velocity uncertainty) in the central 15″, but afterwards increases rapidly. PGC 047273 is a non-BCG, and the only one with regular kinematics (see Section 4.1.2), with a significantly varying δk₀, which also shows a jump from the central velocity within 2″. This galaxy also has a nuclear dust disc, unusually located in a near polar orientation, with LINER ionisation properties (Pagotto et al. 2021). Finally, PGC 047752, PGC 019085 and PGC 099522 are also worth highlighting, with Δk₀ > 20 km/s beyond the effective radius. PGC 099522 is a non-BCG, but its δk₀ is obtained from only half of the field as its velocity map is obstructed by two satellite galaxies (see Appendix A). The extraction of KINEMETRY parameters are therefore less secure and δk₀ is uncertain, and we do not consider it as an k₀ outlier. In summary, two galaxies have a central variation of k₀ (PGC 047273 and PGC 073000), and four galaxies have variable k₀ beyond the effective radii, where PGC 047202 shows the strongest variation from about 0.5 × R_e outwards.

Separating the ATLAS^3D sample into slow and fast rotators shows a potential trend where large variations in the systemic velocity are somewhat more likely for slow rotators (11% slow rotators compared to 7% fast rotators). In the M3G sample, excluding the central k₀ variation, all galaxies with large Δk₀ are BCGs and slow rotators (see Section 4.2).

In Appendix B, we describe three methods to correct M3G velocity maps for the variation of the k₀ term and determine that this variation can be linked with the existence of an additional stellar component. Once the spectral contribution of this components is removed, the new δk₀ is essentially constant. Therefore, we conclude that being massive, a BCG, and a slow rotator favours conditions for existence of multiple kinematics components in the outer regions, likely not in equilibrium with the host galaxy. We discuss this issue further in Section 5.4. Corrected velocity maps are shown with the rest of the kinematics in Fig. C.1.

4.1.2. Non-regularity of M3G velocity maps

After removing the spatially varying k₀ as described in Appendix B, we run KINEMETRY on the velocity maps presented in Fig. C.1. We verified that there are no significant differences between KINEMETRY results obtained on the original or “cleaned” velocity maps, as one would expect due to the orthogonality of the harmonic coefficients. Minor differences are occasionally found in degenerate parameters that define the shape and orientation of the best-fitting ellipses, but they are fully within the uncertainties.

All KINEMETRY parameters were left to vary freely, except the centre, which was fixed to the peak of the surface brightness. In several cases, especially those with very complex velocity maps, the flattening parameter was found not to be well constrained and oscillate over the field, which happens when the assumption on the cosine law is not justified. To constrain the analysis we restricted the ellipse flattening to be equal or larger than the global flattening of the stellar distribution, Q_min ≥ Q_p. Example plots of kinemetric parameters (PA_kin, Q_kin, k₁ and k₁/k₅) for all galaxies are presented in Fig. S3 of the supplementary material. Relevant values are presented in Table S2.

The primary use of KINEMETRY is to quantify irregularities on velocity maps and to highlight kinematic features. We follow the definitions from Krajnović et al. (2011), where the regular rotators (RRs) are galaxies with $\overline{k_5/k_1}\le 0.04$, while non-regular rotators (NRRs) have $\overline{k_5/k_1}>0.04$. The average values are weighted by luminosity and obtained using eq. (1) from Krajnović et al. (2008), see also Ryden et al. (1999).

We applied the above-mentioned classification into kinematic groups a) to e), with a small modification for galaxies that exhibit multiple velocity reversals; they are flagged as multi-spin (MS)¹⁵. Therefore, a galaxy showing two spin reversals is of type MS-2 (e.g. a counter-rotating KDC such as IC 1459, Group c, Franx & Illingworth 1988). Furthermore, when the difference between the PA_kin and the photometric major axis (PA_phot) is close to 90°, the rotation is ‘around’ the long-axis. If such a galaxy has two spin reversals, where one components is in a long-axis rotation, it will be classified as MS-2L, ‘L’ standing for long-axis rotation (e.g. the classic KDC NGC 4365, Davies et al. 2001). For distinction, we refer to the much more common rotation around the minor axis as ‘short-axis rotation’. This nomenclature is based on the connection of the orbital families and the type of rotation they generate: The short-axis tubes rotate around the short (minor) axis of a galaxy, while the long-axis tubes rotate around the long (major) axis of a galaxy¹⁶.

Velocity maps of eight M3G galaxies have two spin reversals (KDCs), and there are at least four galaxies with evidence for more than two kinematic components oriented in different directions. PGC 046832 (MS-5L) is perhaps the most complex example, and we show it in Fig. 5. Looking along the major axis of the galaxy (traced by the solid line on the rightmost panel), one can see three spin reversals, the first within central 3″, the second between ∼5 − 10″ and the third beyond 10″. While these flips between the prograde and retrograde motions are seen in projects, we nominally associate them with to the short-axis rotation. Only detailed modelling (as done in den Brok et al. 2021) can associate visible structures with internal orbital families. Looking along the minor axis (traced by the dashed line) there are two spin reversals, the first ∼3 − 8″ and the second beyond 10″. We associate these with the long-axis rotation reversals. While the orientation of each spin direction is not exactly aligned to the projected major and minor axes, one can count a total of five distinct spin reversals.

Fig. 5. Example of a multi-spin galaxy: PGC 046832. This galaxy has five visible changes of the rotational orientation. Left: Kinematic position angle as a function of radius (1 kpc ∼ 1″). The photometric major axis is shown with a dashed black line, and the median position angles of five components are highlighted with short orange lines. Middle: Observed velocity map. Right: Reconstructed velocity map using low-order KINEMETRY coefficients highlighting only the rotational component of the velocity map. The solid and dashed red lines are the project major and minor axes of the galaxy, respectively. Along the major axis in projection, there are three reversals of the short-axis rotation (prograde, retrograde, prograde), while along the minor axis in projection, there are two (prograde and retrograde), making this galaxy a case with five multi-spin components.

Thirteen galaxies, or about half of the M3G sample, have only one kinematic component (i.e. MS-1). Two of these galaxies have long-axis rotation over the full map, while five galaxies can be classified as RR, with short-axis rotations and nearly constant PA_kin. The velocity maps of other MS-1 galaxies are characterised with kinematic twists (KT) with a steady progression of the PA_kin (both converging to and diverging from the photometric major axis). Long-axis rotation is often found in the innermost component (e.g. PGC 047752, MS-4L), but other options are possible: PGC 047590 (MS-3L) has two prolate-like components that are also counter-rotating with respect to each other, or PGC 007748 (MS-2L) with an outer long-axis rotation. The main characteristic of galaxies with multiple spin components is that all, with a possible exclusion of PGC 046785, have at least one component in (approximate) long-axis rotation. The central component of PGC 046785 is misaligned by about 50^o from the photometric major axis, sufficiently distinct from rotation around the major axis, but strongly misaligned.

The existence of RR galaxies in the M3G sample is quite striking, considering that they are all very massive galaxies. They are, however, only found in non-BCGs, typically fainter galaxies in our sample. In Abell 3558, the brightest RR galaxy (PGC 047177) is ranked as the 4th brightest galaxy in the cluster. In Abell 3556, the second brightest galaxy (PGC 046785) is marginally a RR, with $\overline{k_5/k_1}=0.037\pm 0.004$ within the effective radius, but becomes a NRR at larger scales. There are no RR in Abell 3562, where the luminosity differences between the BCG and other three members are the smallest.

In summary, massive galaxies often have complex stellar kinematics with multiple spin reversals, while the long-axis rotation is common. Nevertheless, some exhibit regular velocity maps consistent with disc-like rotation, but these are typically less bright (massive) galaxies.

4.2. Stellar angular momentum

We follow Emsellem et al. (2011) for calculating the global stellar angular momentum λ_Re within an elliptical aperture of ellipticity ϵ_Re, which has the same area as the circle of radius R_e. The uncertainties were estimated by perturbing the velocity and velocity dispersion maps via Monte Carlo simulations. Figure 6 shows the λ_Re − ϵ_Re diagram for the M3G galaxies, plotted together with the results from the ATLAS^3D (Emsellem et al. 2011) and MASSIVE surveys (Veale et al. 2017b), as well as a selection of nearby brightest group galaxies observed with MUSE (Loubser et al. 2022). We also add the satellites of M3G galaxies for completeness, which are all, as expected, fast rotators. λ_Re values of M3G galaxies are listed in Table 1 and Table S1 of the supplementary material.

Fig. 6. Comparison of λRe versus ellipticity ϵRe for the M3G sample and galaxies drawn from several reference surveys. The M3G galaxies are shown as large red (BCGs) and blue (non-BCGs) circles, with uncertainties smaller than the symbols. Small satellites of M3G galaxies found in the MUSE FoV are represented as teal circles, while small light green circles are for ATLAS3D galaxies. Triangles are used for galaxies from the MASSIVE survey, where filled (brown) triangles are their BGGs, taken from Veale et al. (2017b). Dark red hexagons are BGGs from Loubser et al. (2022). A solid dark green line separates fast rotators from slow rotators, as $0.31\times \sqrt{{ϵ}_{\mathrm{Re}}}$ (Emsellem et al. 2011), while the solid grey line, given by Cappellari (2016) as 0.08 + ϵ/4, is an alternative separatrix. The dashed magenta line shows the edge-on view for ellipsoidal galaxies integrated up to infinity with β = 0.7 × ϵ (Cappellari et al. 2007). The dotted black lines show the location of galaxies, originally on the magenta line, as the inclination is varied by 10°, while the dashed black lines trace locations for galaxies with intrinsic ellipticities from 0.25 (lowest) to 0.85 (highest) in step of 0.1.

There are six fast rotators in the M3G sample, all classified as non-BCGs and RR, which are, as mentioned earlier, lower luminosity galaxies in Abell 3556 and Abell 3558. All M3G BCGs and the remaining five non-BCGs (45%) are classified as slow rotators, which means 76% of all M3G galaxies are slow rotators. Within the M3G sample, Abell 3562 in the SSC has the largest number of slow rotators (four). Notably, slow rotators are found in closer proximity to other more massive galaxies. We show this in Fig. S1 of the supplementary material where we plot the M3G SSC subsample over the distribution of galaxies with spectroscopic redshifts within the SSC (z = 0.048 ± 0.004) from the Haines et al. (2018) catalogue.

Notably, PGC 046785, marginally a RR within the effective radius, has also the lowest λ_Re value among the fast rotators. It remains a fast rotator using the alternative Cappellari (2016) fast versus slow rotator separation. PGC 046785 has a strong kinematic twist of 60° and has three visible multi-spin components. Similarly as it shows significant departures from regular rotation beyond the half-light radius, estimating its λ_R and ϵ values at (2 × R_e) would make it a slow rotator. In contrast, PGC 047355 has the largest λ_Re value in the sample, placing it among one of the fastest rotating ETGs. We also note that PGC 047202, the brightest (most massive) and the largest galaxy in the sample, found in the centre of the richest cluster, also has the largest λ_Re among M3G slow rotators. Finally, PGC 019085 and PGC 048896 are the flattest galaxies in the sample and occupy the sparsely populated region of λ_Re − ϵ_Re with ϵ > 0.4 (Graham et al. 2018).

In Fig. 7 we show the dependence of λ_Re on the absolute brightnesses of galaxies. We also use the calibration from Eq. (2) of Cappellari (2013) to transform the K_s band absolute magnitudes into masses, putting all three samples onto the same mass system. Combining the data from the ATLAS^3D and other surveys of massive galaxies (Veale et al. 2017b; Loubser et al. 2022), one can see a relatively abrupt disappearance of fast rotators above a threshold of M_K ≈ −26 or M ≈ 8 × 10¹¹ M_⊙, where there are essentially no fast rotators (exceptions are NGC 3158 and PGC 046785).

Fig. 7. Distribution of the normalised ${\lambda }_{\mathit{Re}}/\sqrt{ϵ}$ spin parameter as a function of the Ks band absolute magnitude. The horizontal dashed line (at 0.31) separates slow rotators (below) from fast rotators (above). In the plot are galaxies from the M3G (large circles), MASSIVE (triangles), and ATLAS3D (small green circles) samples. The M3G galaxies are split into BCGs (orange) and non-BCGs (blue), while other galaxies are from the same samples as on Fig. 6. The upper axis is in units of log solar masses, obtained with Eq. (2) from Cappellari (2013).

4.3. Correlation with the h₃ Gauss-Hermite moment

In order to investigate what sustains the relatively high λ_R values in some massive galaxies we make use of the Gauss-Hermite h₃ moments, which measure skewness or asymmetric deviations of the LOSVD from a Gaussian distribution. Skewed LOSVDs can be constructed by superimposing two kinematic components with Gaussian LOSVDs, but of different mean values and dispersions. Scorza & Bender (1995) showed that a skewed LOSVD, with a large h₃ value, is typical of a rotating disc. The relation between h₃ and the normalised mean motions, V/σ, is particularly useful, as fast rotators display a strong anti-correlation in a V/σ − h₃ diagram, which is absent for slow rotators (Bender et al. 1994; Krajnović et al. 2011; van de Sande et al. 2017).

Frigo et al. (2019) introduced a quantitative measure for the V/σ − h₃ (anti-)correlation using the following expression:

$$\begin{array}{c}\hfill {\xi }_3=\frac{{\sum }_i{F}_i{h}_{3,i}(V_i/{\sigma }_i)}{{\sum }_i{F}_i{h}_{3,i}^2},\end{array}$$(3)

where the local flux F_i, the mean velocity V_i, the velocity dispersion σ_i and the skewness parameter h_3, i, measured for each spatial bin i within the kinematic maps, provide an estimate of the global slope of the V/σ − h₃ relation. The ξ₃ parameter is tuned such that when V/σ and h₃ are (anti-)correlated, h₃ = (1/ξ₃)V/σ and ξ₃ is the inverse of the slope of the (anti-)correlation. Based on numerical tests (Frigo et al. 2019), fast-rotating galaxies have ξ₃ < −4, while slow rotators are expected to have values close to 0, or mildly positive values. We calculated the ξ₃ values and list them in Table S2 of the supplementary material.

The left panel of Fig. 8 shows the distribution of V/σ − h₃ values for all individual bins which make up the kinematic maps of M3G galaxies. They are divided into those that belong to RR and NRR galaxies (Table 1). The anti-correlation is evident for RR galaxies, while it is barely existent for NRR galaxies, as emphasised by blue dashed and solid orange lines. We also make comparisons between ξ₃ slopes calculated by grouping galaxies in different distributions, for RR and NRR, BCGs and non-BCGs, and based on the distribution of kinematic misalignments in the M3G sample (Krajnović et al. 2018, the classification is repeated in Table S2 for completeness). We also add the ξ₃ values for the RR and NRR galaxies from the ATLAS^3D sample of local ETGs (Krajnović et al. 2011). The comparison indicates that M3G RR, “oblate” (aligned) galaxies, and non-BCGs have more negative values of ξ₃, similar to the value given by the local RR ETGs. M3G NRR, BCGs and “prolate” (highly misaligned) galaxies have ξ₃ values similar to those of local NRR ETGs. The ξ₃ values relating to the slopes of various samples in Fig. 8 are listed in Table S3 in the supplementary material.

Fig. 8. Observed anti-correlation between V/σ and h3. Left: Local V/σ − h3 relation for every spectrum within the effective radius and with an error on h3 smaller than 0.05 for M3G galaxies. Red contours are for NRR galaxies, and blue contours are for RR galaxies, and they are based on logarithmic number counts starting from 0.25 with a step of 0.5. Straight lines show the slopes of V/σ − h3 distributions estimated by the ξ3 parameter (see Section 4.3 for details) for M3G galaxies separated into RR (blue) and NRR (red), BCGs (black solid) and non-BCGs (black dashed), and three classes based on the kinematic misalignment angle (see Table S2): ‘oblate’ (violet dashed), ‘triaxial’ (violet dotted) and ‘prolate’ (violet solid). The values for the local ETGs from the ATLAS3D survey, also separated into RR and NRR galaxies, are plotted in green. Right: Distribution of individual ξ3 values as a function of the normalised spin parameter where the value of 0.31 separates slow rotators from fast rotators (vertical dashed line). The M3G galaxies are divided into BCGs (red) and non-BCGs (blue) and shown with large red and blue symbols, respectively, while the small light green symbols are for local ETGs from the ATLAS3D survey. Dark red hexagons are from Loubser et al. (2022). The error bars on M3G values are typically smaller than the symbols size.

The comparison of the ξ₃ values for the M3G Sample (Table S2), as well as for the ATLAS^3D and the Loubser et al. (2022) samples, is shown in the left panel of Fig. 8. Errors were estimated using Monte Carlo simulations for each kinematic map (V, σ, h₃). As expected, BCGs have low ξ₃, with one galaxy showing a positive value. Non-BCGs also often have low ξ₃ values, including one that was classified as RR and fast rotator (PGC 047177). This galaxy indeed does not have a coherent h₃ map, in spite of a very regular velocity map and a generally high amplitude of rotation (Fig. C.1). Other fast rotators in the M3G sample, however, show a spread of ξ₃ values from ∼ − 2 to −11.

Analysis of the numerical simulations showed that the value of ξ₃ depends on the orbital compositions of a galaxy (Hoffman et al. 2009; Frigo et al. 2019). Strongly negative values of ξ₃ are characteristic for orbital distributions dominated by short-axis tubes, while values close to zero, or even positive, are found in systems with high fractions of long-axis tubes. This is important as it is not only the existence of stellar discs, but the general existence of coherent and dominant streaming around the short axis that results in large negative ξ₃ values. Furthermore, it suggests that the M3G galaxies with high angular momentum contain large fraction of short-axis tubes, which in some cases are completely dominating the orbital distribution (e.g. for ξ₃  − 8.4 the fraction of short-axis tubes is 60%, Frigo et al. 2019).

5. Discussion

The main finding of this work is a surprising richness of features on velocity maps of galaxies with masses above ∼10¹² M_⊙. High quality IFU data reveal velocity maps with multiple-spin reversals, not seen before. On the other hand, regular rotation and fast rotators persist even among very massive galaxies.

5.1. Multi-spin nature of velocity maps of massive galaxies

The M3G sample comprises galaxies that do not show net rotation (Group a: 2/25), galaxies with complex irregular velocities (Group b: 9/25), galaxies with multiple spin reversals (Group c: 8/25), and galaxies with regular rotation (Group e: 6/25). Noteworthy is that 44% of Group b and all galaxies of Group c have at least one component in long-axis rotation (52% of the full sample). Recent simulation studies predict that the emergence of prolate morphology and long-axis rotation is linked to the formation of massive galaxies (Ebrová & Łokas 2017; Li et al. 2018; Lagos et al. 2022). M3G data contribute by showing that the long-axis rotation is prevalent among massive galaxies, in the sense that it exists even if it is not a dominant component.

The velocity maps of the M3G Group c galaxies are remarkable for atypically large number of kinematic features. Galaxies up to 5 − 6 × 10¹¹ M_⊙ generally exhibit at most two different stellar kinematic components, with respective spins aligned with either the principal major or minor photometric axis (e.g. NGC 4365, the quintessential example of a KDC galaxy; Davies et al. 2001). While KDC galaxies are also present in the M3G sample of galaxies, we reveal here systems that have much more intricate velocity maps, as well as multiple spin reversals. The most remarkable cases are the following:

• PGC 019085 (MS-3L), with a central KDC, outer counter-rotation, and an additional component with long-axis rotation;

• PGC 046832 (MS-5L), presented earlier (Fig. 5) with five spin reversals;

• PGC 047590 (MS-2L), with prograde and retrograde long-axis rotations; and

• PGC 047752 (MS-4L), with a central 5″ region showing both short and long-axis rotations and a retrograde short-axis rotation at larger radii, which twists strongly but incompletely towards long-axis rotation at the edge of the MUSE field.

We show these systems (except PGC 046832) in Fig. 9, illustrated by their observed velocity maps and KINEMETRY-reconstructed velocity maps (see also Figs. C.1 and S3 for more details). Another potential case is PGC 043900 (see KINEMETRY-reconstructed velocity map in Fig. S3). It is classified as a non-rotating Group a galaxy. Despite the low speeds, when using the harmonic analysis and averaging the velocities along the ellipses, one can see hints of both prograde and retrograde short-axis rotations throughout the main body of the galaxy, as well as evidence for long-axis rotation confined to the central 10″.

Fig. 9. Galaxies from M3G with extreme velocity maps. The observed map is on the left, and the right panel shows the reconstructed velocity maps using low-order KINEMETRY coefficients, highlighting only the rotational component. We note that the colour bars have somewhat different scales and have been tuned to highlight the features on individual maps. Solid and dashed red lines are the projected major and minor axes of the galaxy, respectively. Kinematic features are best seen by following the spin reversals along the major and minor axes. Another extreme, PGC 046832, was already shown in Fig. 5, which has five visible changes of the rotational orientation.

Complex velocity maps of equilibrium models are possible only in non-axisymmetric objects, as the observed velocity map is a luminosity-weighted average of the stellar orbits that build the galaxy (Statler 1987). In addition, the shape of the velocity map also depends on the projection, or the viewing angle, with respect to the orientation of the galaxy (Franx 1988; Statler 1991; Arnold et al. 1994). In his analysis of triaxial models based on the “perfect ellipsoid” – a Stäckel potential form (de Zeeuw 1985) – Statler (1991) demonstrated that strong kinematic twists, long-axis rotation, and KDCs arise naturally from the presence of three distinct families of tube orbits and the effects of varying viewing angles. This is true even when rotation is restricted to only prograde streaming along the short- and long-axis tubes. Due to projections effects, the Statler (1991, Fig. 3) velocity maps included classic KDCs of the NGC 4365 type, as well as predicted the possibility of observing types of long-axis counter rotation present in PGC 047590 (Fig 9). The velocity maps of real galaxies are, however, more complex than the ‘perfect ellipsoid’ models with streaming around two axes, because we observe both prograde and retrograde rotations.

Evidence for stellar counter-rotation in oblate disc galaxies is abundant (starting from the archetypal NGC 4550; Rubin et al. 1992). They can be traced to the accretion of gas on retrograde orbits (e.g. Bertola et al. 1992; Merrifield & Kuijken 1994; Pizzella et al. 2004; Coccato et al. 2013), mergers (Bois et al. 2011), or potentially to cold streams (Algorry et al. 2014). In massive ETGs, which rotate more slowly, it is much more difficult to directly extract the kinematics of different components by spectroscopically disentangling them. However, one can investigate the internal structure of ETGs using dynamical models, especially with orbit-based modelling (Schwarzschild 1979). A conceptual proof of this approach was presented with a dynamical model of NGC 4550 (Cappellari et al. 2006), demonstrating that the observed kinematics can be reproduced by placing ∼50% of the mass into the counter-rotating component. Another case was the classic KDC in NGC 4365, which features two orthogonal spin axis coinciding with the major and minor photometric axes. Using triaxial dynamical models, van den Bosch et al. (2008) revealed that to reproduce the kinematics, the internal structure requires three different orbital components with streaming: prograde and retrograde short axis tubes and (prograde) long axis tubes¹⁷. Further studies using Schwarzschild’s method showed that galaxies with KDCs or counter-rotating discs can be similarly reproduced by superposition of orbits (e.g. IC1459, NGC 4473, NGC 5813 by Cappellari et al. 2002, 2007; Krajnović et al. 2015; Mulcahey et al. 2021). These studies showed that even what is observed as non-rotation can be explained as a (luminosity or mass weighted) combination of prograde and retrograde orbits.

One of the most complex-looking M3G galaxies, PGC 046832 (Fig. 5), was modelled in a similar way with triaxial Schwarzschild models by den Brok et al. (2021), who concluded this galaxy is built of prograde and retrograde short axis tubes (42% of the mass), prograde and retrograde long axis tubes (34%) and box orbits (23%). Crucially, all four families of streaming orbits are needed to reconstruct the five spin reversals observed in the velocity map. Therefore, it is reasonable to assume that all other multi-spin galaxies will require a similar orbital set-up to reproduce their velocity maps, in addition to the projection effects due to their triaxial nature and specific viewing angles. A working hypothesis is that the orbital families permeate through each other and the full body of the galaxy, while the features in velocity maps depend on the relative spatial contributions of various orbital families, their prograde and retrograde fractions, the viewing angles and projection effects (Statler 1991; Arnold et al. 1994). It follows that features visible in the velocity maps (e.g. a KDC) are not a specific well defined entities in real space, in the sense of being a physical and spatially confined component, separated from the rest of the galaxy body (van den Bosch et al. 2008; Krajnović et al. 2015; Mulcahey et al. 2021).

The simplest explanation of retrograde motions (counter-rotation) is that it originates from accretion of external material. In case of oblate counter-rotating discs, the accretion of gas is sufficient (Algorry et al. 2014). For triaxial and massive galaxies, strong tidal interactions and mergers are the most likely explanation for the observed richness of kinematics (Hoffman et al. 2010; Ebrová et al. 2021), especially when they are gas-free mergers of similar mass galaxies (Röttgers et al. 2014; Naab et al. 2014; Frigo et al. 2019). This is supported by the variety and apparent randomness of the observed velocity maps: velocity maps of massive galaxies rarely exhibit significant similarities, even when the galaxies have the same number of spin reversals. Furthermore, multi-spin velocity maps, and especially those with more than two spin reversals, occur in brighter (more massive) galaxies, which are more likely to experience gas-free mergers. Finally, the fine-tuned balance among orbital families that might even lead to the cancellation of detectable streaming motions requires several merging events. Such systems are the rarest (3% in ATLAS^3D survey and < 1% in M3G), and happen in the most massive galaxies, which again are more likely to experience multiple minor and major mergers (Naab et al. 2014).

Conversely, velocity maps of RR galaxies resemble each other to a much higher degree, which makes sense if their formation is predominantly through dissipative processes, which include star formation or and disc formation. Their subsequent passive evolution might transform them into massive fast rotator ETGs, but not into multi-spin slow rotators. Therefore, multi-spin structures on velocity maps are possible only if galaxies are triaxial, allow for streaming along the long and the short axis, and are formed through interactions and mergers.

5.2. Slow rotation at the highest masses

All M3G BCG are slow rotators. Additionally, all non-BCG M3G members of Abell 3562 cluster are slow rotators, as well as the two brightest non-BCGs in Abell 3558. Abell 3562 is somewhat special as its four brightest members have similar brightness (e.g. there is no typical brightness gap between the BCGs and non-BCGs, Lauer et al. 2014), the cluster does not seem settled and several galaxies are bright X-ray sources (see Section 2.1). Still, the M3G sample follows a trend: slow rotators are BCGs, possibly include one or two highest ranked galaxies in a cluster, and they are found in regions more densely populated with other galaxies (in comparison with fast rotators). This is consistent with previous studies (Cappellari 2013; Jimmy et al. 2013; Fogarty et al. 2015; Oliva-Altamirano et al. 2017; Brough et al. 2017; Graham et al. 2019b), where BCGs of massive clusters are exclusively slow rotators, and only a few other galaxies within the cluster, typically at centres of sub-clusters groups, can also be slow rotators (D’Eugenio et al. 2013; Fogarty et al. 2014; Scott et al. 2014; Graham et al. 2019a).

In this respect, our results support the existence of the kinematic morphology - density relation (Cappellari et al. 2011b; Cappellari 2013). While the environment might not be the main driver (Brough et al. 2017; Greene et al. 2017; van de Sande et al. 2021a), the mass function of slow rotators depends on the environment (Graham et al. 2019b). The M3G sample contains some of the largest and brightest galaxies and extends the current sample of massive galaxies to K_s ≈ −27 absolute magnitude. As a consequence we report that galaxies brighter than −26 mag are almost exclusively slow rotators (Fig. 7). This level of brightness is only achieved by BCGs, which also do not belong to the general population (or the luminosity function) of massive galaxies (e.g. Lin et al. 2010). The M3G sample illustrates the interconnectedness of galaxy mass growth and the loss of angular momentum. Slow rotation of BCGs is naturally explained through repeated gas-poor major mergers, as well as the frequent accretion of smaller satellites at their “privileged” position at the bottom of the cluster potential well. Such mergers contribute to the reduction of angular momentum, growth of galaxy mass and diversification of the types of orbits, and explain the variety of kinematic features on velocity maps.

5.3. Fast rotation at the highest masses

Using KINEMETRY we established that there are six galaxies which are consistent with regular rotation. All these galaxies can unambiguously be classified as fast rotators. Their existence is not entirely surprising: among the nearby ETGs, regular rotation is found in galaxies as massive as 5 − 6 × 10¹¹ M_⊙ (Krajnović et al. 2020) and the MASSIVE survey reports fast rotators up to 10¹² M_⊙ (Veale et al. 2017a). None of the M3G fast rotators is a BCG. Instead, they are by brightness all lower ranked galaxies in their clusters. The only exception is the second ranked galaxy in Abell 3556, PGC 046785, which, as argued before, is marginally a fast rotator within its effective radius. All other M3G fast rotators have regular velocity maps, suggesting the existence of a single dominant orbital family of short-axis tubes.

This assessment is based on ξ₃ values of M3G fast rotators, which are all, except one, strongly negative, implying an anti-correlation V/σ − h₃ (Fig. 8, left). Simulations predict that such galaxies have high fraction of short-axis tubes, dominating the orbital distribution, and are likely to have experienced gas rich mergers producing fast rotators (Hoffman et al. 2009; Naab et al. 2014; Röttgers et al. 2014; Frigo et al. 2019). An exception is PGC 047177, which therefore must have a more complex internal orbital structure, in spite of a regular velocity map.

Fast rotators are galaxies less bright than −26 mag (K_s band), and approximately less massive than 8 × 10¹¹AR > M_⊙. As there is a clear relation between the formation of slow rotators and mergers of certain type (Hoffman et al. 2010; Bois et al. 2011; Naab et al. 2014), this brightness (mass) limit possibly defines an approximate threshold above which formation of fast rotators is less likely.

5.4. BCGs and their massive companions

There is ample evidence that BCGs are different from other galaxies (e.g. Tremaine & Richstone 1977; Graham et al. 1996; Lauer et al. 2014; Loubser et al. 2018, 2020), but there are also similarities between BCGs and the highest mass non-BCGs. Kinematically, for example, short-axis rotation is found both in BCGs and non-BCGs, and there are also BCGs that have very small kinematic misalignment angles (Krajnović et al. 2018; Ene et al. 2018). We showed that there are non-BCGs that are slow rotators and have multiple-spin reversals, but they are rare. Regular rotation and fast rotators are only found in non-BCGs, and long-axis rotation predominantly occurs in BCGs: 10/14 (71%) compared to 3/11 (27%) in non-BCGs.

We also presented another difference between BCGs and non-BCGs: a radially varying k₀ term. As discussed in Section 4.1.1 and Appendix B, in most cases this is likely an indication of the existence of another kinematic component at outer radii. We do not present the evidence here and postpone doing so for a future paper. Galaxies that show a large δk₀ variation also have radially raising velocity dispersion profiles (see e.g. velocity dispersion maps of PGC 015524 and PGC 047202 on Fig. C.1). Rising velocity dispersion profiles are found in massive galaxies (Veale et al. 2018) and are specific for BCGs and cD galaxies (Brough et al. 2007; Newman et al. 2013) where velocity dispersion reaches the value of the velocity dispersion of the cluster galaxies. An example is NGC 6616 showing radial variation in the mean velocity and the velocity dispersion, both of which approach the mean cluster values (Bender et al. 2015). The authors explain such high velocity dispersion values by assuming that outer stellar component of the BCGs is build up from stars accreted in minor mergers or stripped from cluster galaxies. Our data do not extend to large enough radii, but we see a similar trend where the velocities, parametrised by V_sys + k₀, rise towards the mean cluster velocity values. Therefore, it is possible that the δk₀ variations trace an increasing contribution of a component made up of externally accreted stars. As this component likely originates from tidally stripped material of other cluster galaxies, it is (nearly) at rest with respect to the cluster mean velocity and velocity dispersion based on the motions of all cluster galaxies. In this sense, the kinematics of this components could be different from the main body of the BCG, which has its specific location and motion within the cluster.

Numerical simulations (Han et al. 2024) find the increase in the velocity dispersion at outer radii to be related to the high fraction of stars formed ex situ and higher velocity anisotropy which these stars carry. The same authors argue that the variety of the outer velocity dispersion profiles in galaxies reflects their diverse assembly histories. We will investigate the velocity dispersion radial profiles and their connection to δk₀ in an upcoming study. Nevertheless, we conclude that if one can establish a connection between the radial k₀ increase and accretion of stars born elsewhere, which form an additional kinematic component, this would provide further evidence of different formation paths for BCGs and non-BCGs.

6. Conclusions

The observed stellar kinematics are a reflection of the richness of the underlying orbital distributions that build up galaxies. This can be seen in the orderly velocity maps of RRs, which comprise the majority of galaxies, and in the kinematic complexity observed in the velocity maps of the most massive systems. This interpretation relies on the hypothesis that galaxies are built of several orbital families, such as the short-axis and long-axis tubes, rotating around the minor and the major axis of the galaxy, respectively, and that these families leave their imprint on the velocity maps. There are also orbits that have no angular momentum, and although they do not leave an imprint on the velocity maps, they participate in shaping the even kinematical moments, such as galaxy morphology and the velocity dispersion.

In this work we have described the observational setup and data associated with the M3G Survey, a MUSE GTO program, targeting several of the brightest (most massive) galaxies within ∼200 Mpc. From observations of 25 such galaxies, divided into 14 BCGs and 11 non-BCGs, we extracted stellar kinematics and analysed their velocity and h₃ maps. We demonstrated that velocity maps of massive galaxies can show several reversals in the orientation of the visible rotation, notably many more compared to what was known for less massive galaxies (< 5 × 10¹¹ M_⊙), including classical KDCs. In the extreme case of PGC 046832, there are five different spin-flips in projection. Furthermore, the majority of the BCGs and several luminous non-BCGs have at least one component that shows long-axis rotation (rotation around the major photometric axis). These complex velocity maps reinforce the notion that massive galaxies have non-axisymmetric shapes and contain several orbital families, but stars need to stream in both prograde and retrograde directions within their parent orbital families.

We find that there are massive galaxies that have very simple velocity maps (regular rotation) and that they are all fast rotators. They are non-BCGs and also low ranked (in luminosity) within their clusters, in spite of all sample galaxies being among the brightest in the nearby Universe. Fast-rotating galaxies also show a strong anti-correlation between V/σ and h₃ parameters, which in comparison with numerical simulations suggests their orbital distribution is simple and dominated by short-axis tubes only.

Some BCGs show kinematic evidence for the existence of an additional component not fully synchronised with the main component. The difference in the kinematics points to an external origin of this component, possibly being made of stars stripped from other cluster galaxies. This is only a tentative conclusion, and it will be explored more in a future paper discussing the velocity dispersion and the h₄ moment maps of M3G galaxies.

Overall, our findings support the general paradigm of the formation of massive galaxies (e.g. Oser et al. 2010; Cappellari 2016; Naab & Ostriker 2017) and provide further constraints that the formation of the most massive galaxies is governed by merging. The mergers are typically gas-poor, but with diverse characteristics such as for the mass ratios, orientation of the orbital angular momentum, and frequency depending on the local environment. Velocity maps remain a very good indicator of the complexity of evolutionary histories, especially at the tip of the galaxy mass distribution.

Data availability

The supplementary material is accessible at https://doi.org/10.5281/zenodo.18459530

Acknowledgments

We thank Ilaria Pagotto and Mark van den Brok for contributions to data analysis and discussions at early stages of the project. DK acknowledges support from the grant: KR 4548/4-1 of the Deutsche Forschungsgemeinschaft. L.A.B. acknowledges support from the Dutch Research Council (NWO) under grant VI.Veni.242.055 (https://doi.org/10.61686/LAJVP77714). We have used the services from SIMBAD, NED, and ADS as well as the Python packages NumPy (Oliphant 2007), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), and Astropy (Astropy Collaboration 2018), for which we are thankful. The MUSE observations were taken within the following observing programmes: 094.B-0592, 095.B-0127, 096.B-0062, 097.B-0776, 098.B-0240,099.B-0148,099.B-0242, 0102.A-0327.

References

Appendix A: Extraction of satellites and their kinematics

Figure S2 in the Supplementary Material presents three examples of M3G galaxies with fore- or back-ground satellites. Plots for each galaxy are distributed over two columns, where the first row shows the map of the χ² used to allocate whether a bin belongs to the main galaxy, satellite galaxy or both, and the distribution of these bins (see Section 3.1 for details). PGC 047355 is the only case in the sample where both the galaxy and the satellite kinematic components can be disentangled at all bins. In all other cases, the central regions of the satellites are completely dominating its light and we are not able to obtain any information about the main galaxy kinematics. The spectra in which we can separate the components are typically confided to a small region at the boarders of the satellites. An exception is PGC 099188*, used as an example in the main text (Section 3.1), where the multi-component pPXF fit is able to separate both the satellite and the main galaxy component over the most of the field, except in the central 10″ of the satellite. Another extreme is shown by PGC 065588 and its satellite in the lower left corner, where only a handful of bins at the edges of the satellite can be used to extract both kinematic components, while for PGC 099522 the two satellites completely dominate one side of the galaxy and the main galaxy kinematics cannot be extracted. Nevertheless, in most cases, the two component fitting procedure yielded a much increased region with kinematics of the main galaxy. Furthermore, it provided a natural and an optimal way of finding bins at the edges of the satellites that still contain significant fraction of the main galaxy light to extract its kinematics.

There are 13 M3G galaxies with satellites whose kinematics we attempt to separate from the main galaxy. For six M3G galaxies, there is only one satellite in the MUSE FoV. Five M3G galaxies have two interlopers, while PGC 003342 galaxy has three satellites, and PGC 047202 has seven. PGC 047202 is the largest galaxy in the sample and found in the centre of the densest cluster, but all of the seven satellites that we separate kinematically are found outside of the central 1 arcmin (central MUSE pointing).¹⁸ Therefore, the reason why we detect such a larger number of satellites compared to the rest of the sample is due to the fact that we mosaic this galaxy with a 2 × 2 MUSE pointings quadrupling the area compared to other galaxies. The location of satellites can be seen on Fig. C.1. In total, we extracted kinematics for 26 satellites, of which half was in the foreground and half in the background of the main galaxies (when averaged over all M3G sample galaxies). We do not show the kinematics of satellites (see Bacon et al. 2023, for kinematics of PGC 047202 satellites), but their kinematics is dominated by rotation, and they are all classified as fast rotators (Fig. 6). In Table S1 of the Supplementary Material we provide an overview of satellites and their location on the MUSE maps. We also provide their ellipticities and λ_Re values shown in Fig. 6, obtained using the same method as for the main galaxies (Section 4.2).

Appendix B: Correcting velocity maps for the effect of variable k₀ parameter

The velocity measured by the pPXF is a combination of the peculiar velocity, V_pec(x,y), specific for each bin and the velocity of the barycentre of the galaxy. The former is the velocity that describes the internal dynamics of the galaxy, and is directly related to the orbital distribution, while latter is a measure of the cosmological redshift of the galaxy (see eq. (3) and Section 2.2 in Cappellari 2023). This is usually referred to as the systemic velocity, V_sys. For clarity, we use an example of a two dimensional velocity map of a disc galaxy, described as

$$\begin{array}{c}\hfill V_{\mathrm{pec}}(x,y)=V(R,\theta )=V_{\mathrm{sys}}(R)+V_{\mathrm{rot}}cos(\theta )sin(i),\end{array}$$(B.1)

where i is the inclination of the disc with respect to the plane of the sky (for i = 90° the galaxy is “edge-on"), θ is the azimuthal angle measured from the major axis of length R (θ = 0°), and V_rot is the amplitude of the rotation. The total amplitude of V_rot depends on the galaxy mass (i.e. traces circular velocity in a disc galaxy) and is modulated by the viewing angle, while its change across the field of the galaxy follows a cosine variation. If the galaxy redshift is precisely known, V_sys = 0, but this is rarely the case, and one measures V_sys ≠ 0 and a certain error ΔV_sys. The systemic velocity is not related to the internal kinematics or dynamics, and it should be constant across the galaxy. However, variations of measured V_sys for each V_pec are expected, and in practice one estimates a global value of V_sys as a median of the velocities across all independent spatial bins.

The variation of k₀ across the field, as defined in KINEMETRY is, in principle, related to the uncertainties of the V_sys, and its average value also provides an estimate for V_sys. However, if the kinematics is not regular (a case different from the ideal disc of eq. B.1), additional contribution might be present, tracing anther kinematic component. This component could potentially not be fully bound to the main galaxy (e.g. following a recent merger or accretion), and retains a “memory" of the cluster kinematics (Bender et al. 2015). In the outer regions where the main galaxy components becomes fainter, the contribution of the secondary components can bias the extracted kinematics, potentially detectable as a systematic shift in mean velocity. A natural way to detect this would be through a variation of k₀.

Here we first present a way to remove the variation of k₀ across the field and determine the underlaying velocity of the main galaxy component, and subsequently we perform two different pPXF extractions to test the hypothesis that the variation of k₀ is due to existence of another kinematic component, detectable at large radii (beyond R_e as suggested by Fig. 4).

In order to determine the amplitude of k₀ variations, we extracted velocities along the fits to the isophotes, and not the best fitting ellipses (for a velocity map) as is usually done. In this way we trace the shape of the stellar distribution, and not the properties of the velocity map. We first run KINEMETRY on the MGE models of M3G galaxies (the same as used for masking in Sections 3). KINEMETRY is run in its photometric mode (all harmonic terms) to determine the photometric position angle (PA_p) and flattening (Q_p). The MGE models are used as smooth representation of galaxies, with slowly varying Q_p and PA_p, without satellite galaxies or compact objects.

The next step was to run KINEMETRY on the velocity maps, imposing the PA_p and Q_p in the extraction, and measuring the value of k₀, via harmonic analysis. We note that the ATLAS^3Dk₀ values, compared in Fig. 4, were extracted in a slightly different way by running KINEMETRY on velocity maps only. The middle panels in Fig. B.1 show the maps of k₀, demonstrating a strong radial dependence, where k₀ rises towards the edge of the MUSE FoV, indicating a strong difference between the inner and outer regions of galaxies. The maps were subsequently subtracted from the original velocities producing the corrected velocity maps. Their examples are presented in the third column of Fig. B.1, and show detailed velocity structures and clear improvements with respect to the original extractions. These maps are used in all subsequent analysis presented in the paper (Fig. C.1).

Fig. B.1. Effects of removing the variable k0. They are best visible in PGC 015524 (top) and PGC 019085 (bottom). Left: Original velocity map extracted using the pPXF from a binned MUSE cube. Middle left: Maps of the k0 term in the harmonic analysis obtained with KINEMETRY on the original velocity maps. Note the different colour scales compared to other plots. Middle right: Corrected velocity map obtained by subtracting the k0 value from each bin from the original velocity map. Right: Velocity maps obtained through pPXF extraction as in the leftmost panels but after first removing a representative halo spectrum (see text for details). Only bins with S/rN> 7 are plotted.

As k₀ maps show the strongest difference at the outskirts of the MUSE field, we tested a hypothesis that this variation is caused by existence of another kinematic component, detectable at large radii. If this secondary component exists, by removing its representative spectrum we should be able to recover the kinematics of the primary component. For both galaxies in Fig. B.1, we selected spectra from 5 − 6 regions, within 3″ apertures, at the very edge of the MUSE field corresponding to high k₀ values (at least one in each quadrant of the field). We fit these outer spectra with pPXF and find their systemic velocities to be substantially different from that of the global galaxy spectrum (extracted within 1R_e). For PGC 015524 the difference reaches 50 km/s and for PGC 019085 about 35 km/s, comparable to the range of the respective k₀ values. For each galaxy, we made a median stack spectrum based on the outer spectra, and subtracted it from each spectrum of the binned data cube. The pPXF fit on the whole data cube was executed as before, now extracting only two moments, V and σ. In the last column of Fig. B.1, we plot the resulting velocity maps for the two galaxies. The outermost bins are masked as the S/rN of the fit was low (< 7). Nevertheless, both maps show clear velocity structures, essentially the same as those obtained by subtracting the k₀ parameter values.

This test suggests that the origin of the large scales variation of k₀ is very likely due to another kinematic component which becomes dominant at large radii. We further tested this hypothesis by setting up pPXF to fit for two kinematics components, of different velocities and velocity dispersions, similar as in the case of fitting satellite galaxies (Section 3.1). We performed this test on PGC 015524 only, as the difference in kinematics between the components are small and large S/N is required for more robust results. After the pPXF extraction, we run KINEMETRY on the resulting velocity maps, estimated k₀ and δk₀ = |k₀(0)−k₀(R)|, and compared with the original δk₀ (e.g. in Fig. 4). We also run KINEMETRY on the velocity map of the primary kinematic components from previous test (map on Fig. B.1), and show all three δk₀ curves on Fig. B.2.

Fig. B.2. Comparison of the absolute radial variation of k0 from three kinematic extractions for PGC 015524. The blue line shows δk0 from the nominal extraction (same as in Fig. 4). The orange dashed line shows δk0 of the main kinematic components from a two components pPXF fit. The dashed dotted green line shows the δk0 from the kinematic extraction after removing a representative spectrum of the secondary component from the cube. The variation in δk0 seems to originate in the existence of a secondary component which dominates the outer parts.

The results of these tests are very suggestive that the variation in k₀ indeed originates in the existence of a secondary component which dominates the outer parts. The secondary components might not be in the equilibrium with the central (main galaxy) component (see discussion in Section 5.4), and could originate from accretion of satellites that distribute their stars at large radii and in the halo of the galaxy.

Appendix C: Kinematic maps

Fig. C.1. Kinematic maps of M3G galaxies. From left to right, the mean velocity, velocity dispersion, h3 and h4 moments are shown. The values (in kilometres per second for the mean velocity and velocity dispersion) in lower-left corners correspond to maximum and minimum of the colour bars on the right of each plot. The red ellipse has the major axis length equal to Re and the median flattening of the galaxy within the effective radius (see Table 1). It shows the elliptical aperture used to extract global effective properties (e.g. σe, λRe, average values of kinemetry parameters, etc).

Fig. C.1. Continued.

Fig. C.1. Continued.

Fig. C.1. Continued.

Fig. C.1. Continued.

All Tables

All Figures

Fig. 1. Properties of the M3G sample. Top left: Histogram of the distances to the BCGs (orange) and the satellites galaxies (blue). Bottom left: Richness of the host galaxy cluster versus the absolute Ks magnitude of the corresponding BCG. BCGs in the SSC are shown as pentagrams (Abell 3558, Abell 3556 and Abell 3562 in order of decreasing brightness of the BCGs). The Virgo and Coma cluster (squares) have been added as references for the ATLAS3D and the MASSIVE Surveys. Right: Size–luminosity relation for the M3G sample galaxies, split into BCGs (red) and satellites (blue). Galaxies from two other surveys have been added for comparison: ETGs from the ATLAS3D survey (light green large squares), spirals from the ATLAS3D survey parent sample (dark green small squares), and ETGs from the MASSIVE survey (brown triangles). ATLAS3D sample is a magnitude-limited sample within 42 Mpc, while the galaxies in the MASSIVE are the brightest galaxies within 100 Mpc. The solid and dashed black lines are relations for fast rotators and S0-Sa galaxies from the ATLAS3D survey, respectively. The red line denotes the so-called ‘zone of avoidance’. All three relations are taken from Cappellari et al. (2011a).

In the text

	Fig. 2. Distribution of the VLT auto-guider seeing measurements averaged over each exposure taken with MUSE for all M3G galaxies. Horizontal lines denote the two seeing limits imposed on the programme. Blue squares show the mean value of all observations for each galaxy. Red circles are exposures targeted to be done in ‘good seeing’ conditions (FWHM < 0.8″).
In the text

Fig. 3. Spectral fits using pPXF for three different spectra from the PGC 099188* MUSE field, which features a large interloping satellite. From top to bottom, the panels show the central spectrum of PGC 099188*, a spectrum from the midpoint between the centre of PGC 099188* and the interloping satellite, and the central spectrum of the interloping satellite. In each panel, the black line shows the observed spectrum, and the red line shows the pPXF fit, with the residual of the fit in green points at the bottom. Green shaded polygons are masked regions excluded from the fit. The middle panel has two additional spectra shown by light and dark blue colours, belonging to the main galaxy and the satellite, respectively. These are model spectra vertically shifted for an arbitrary amount, obtained from the pPXF fits highlighting the satellite and the galaxy contributions to the total spectrum.

In the text

Fig. 4. Radial profiles of the absolute k0 variation for M3G galaxies (δk0 = |k0(0)−k0(R)|) compared to the same from the ATLAS3D galaxies. The M3G galaxies are divided into BCGs (red lines) and non-BCGs (blue lines), while ATLAS3D galaxies are divided into fast (light blue dashed lines) and slow (light red dashed lines) rotators. The regions shaded purple show the standard deviation of all ATLAS3D and M3G (darker shaded regions) velocity errors. Names highlight specific galaxies with the largest difference between the maximum and the minimum values of k0: Δk0 > 20 km/s.

In the text

Fig. 5. Example of a multi-spin galaxy: PGC 046832. This galaxy has five visible changes of the rotational orientation. Left: Kinematic position angle as a function of radius (1 kpc ∼ 1″). The photometric major axis is shown with a dashed black line, and the median position angles of five components are highlighted with short orange lines. Middle: Observed velocity map. Right: Reconstructed velocity map using low-order KINEMETRY coefficients highlighting only the rotational component of the velocity map. The solid and dashed red lines are the project major and minor axes of the galaxy, respectively. Along the major axis in projection, there are three reversals of the short-axis rotation (prograde, retrograde, prograde), while along the minor axis in projection, there are two (prograde and retrograde), making this galaxy a case with five multi-spin components.

In the text

Fig. 6. Comparison of λRe versus ellipticity ϵRe for the M3G sample and galaxies drawn from several reference surveys. The M3G galaxies are shown as large red (BCGs) and blue (non-BCGs) circles, with uncertainties smaller than the symbols. Small satellites of M3G galaxies found in the MUSE FoV are represented as teal circles, while small light green circles are for ATLAS3D galaxies. Triangles are used for galaxies from the MASSIVE survey, where filled (brown) triangles are their BGGs, taken from Veale et al. (2017b). Dark red hexagons are BGGs from Loubser et al. (2022). A solid dark green line separates fast rotators from slow rotators, as $0.31\times \sqrt{{ϵ}_{\mathrm{Re}}}$ (Emsellem et al. 2011), while the solid grey line, given by Cappellari (2016) as 0.08 + ϵ/4, is an alternative separatrix. The dashed magenta line shows the edge-on view for ellipsoidal galaxies integrated up to infinity with β = 0.7 × ϵ (Cappellari et al. 2007). The dotted black lines show the location of galaxies, originally on the magenta line, as the inclination is varied by 10°, while the dashed black lines trace locations for galaxies with intrinsic ellipticities from 0.25 (lowest) to 0.85 (highest) in step of 0.1.

In the text

Fig. 7. Distribution of the normalised ${\lambda }_{\mathit{Re}}/\sqrt{ϵ}$ spin parameter as a function of the Ks band absolute magnitude. The horizontal dashed line (at 0.31) separates slow rotators (below) from fast rotators (above). In the plot are galaxies from the M3G (large circles), MASSIVE (triangles), and ATLAS3D (small green circles) samples. The M3G galaxies are split into BCGs (orange) and non-BCGs (blue), while other galaxies are from the same samples as on Fig. 6. The upper axis is in units of log solar masses, obtained with Eq. (2) from Cappellari (2013).

In the text

Fig. 8. Observed anti-correlation between V/σ and h3. Left: Local V/σ − h3 relation for every spectrum within the effective radius and with an error on h3 smaller than 0.05 for M3G galaxies. Red contours are for NRR galaxies, and blue contours are for RR galaxies, and they are based on logarithmic number counts starting from 0.25 with a step of 0.5. Straight lines show the slopes of V/σ − h3 distributions estimated by the ξ3 parameter (see Section 4.3 for details) for M3G galaxies separated into RR (blue) and NRR (red), BCGs (black solid) and non-BCGs (black dashed), and three classes based on the kinematic misalignment angle (see Table S2): ‘oblate’ (violet dashed), ‘triaxial’ (violet dotted) and ‘prolate’ (violet solid). The values for the local ETGs from the ATLAS3D survey, also separated into RR and NRR galaxies, are plotted in green. Right: Distribution of individual ξ3 values as a function of the normalised spin parameter where the value of 0.31 separates slow rotators from fast rotators (vertical dashed line). The M3G galaxies are divided into BCGs (red) and non-BCGs (blue) and shown with large red and blue symbols, respectively, while the small light green symbols are for local ETGs from the ATLAS3D survey. Dark red hexagons are from Loubser et al. (2022). The error bars on M3G values are typically smaller than the symbols size.

In the text

Fig. 9. Galaxies from M3G with extreme velocity maps. The observed map is on the left, and the right panel shows the reconstructed velocity maps using low-order KINEMETRY coefficients, highlighting only the rotational component. We note that the colour bars have somewhat different scales and have been tuned to highlight the features on individual maps. Solid and dashed red lines are the projected major and minor axes of the galaxy, respectively. Kinematic features are best seen by following the spin reversals along the major and minor axes. Another extreme, PGC 046832, was already shown in Fig. 5, which has five visible changes of the rotational orientation.

In the text

Fig. B.1. Effects of removing the variable k0. They are best visible in PGC 015524 (top) and PGC 019085 (bottom). Left: Original velocity map extracted using the pPXF from a binned MUSE cube. Middle left: Maps of the k0 term in the harmonic analysis obtained with KINEMETRY on the original velocity maps. Note the different colour scales compared to other plots. Middle right: Corrected velocity map obtained by subtracting the k0 value from each bin from the original velocity map. Right: Velocity maps obtained through pPXF extraction as in the leftmost panels but after first removing a representative halo spectrum (see text for details). Only bins with S/rN> 7 are plotted.

In the text

Fig. B.2. Comparison of the absolute radial variation of k0 from three kinematic extractions for PGC 015524. The blue line shows δk0 from the nominal extraction (same as in Fig. 4). The orange dashed line shows δk0 of the main kinematic components from a two components pPXF fit. The dashed dotted green line shows the δk0 from the kinematic extraction after removing a representative spectrum of the secondary component from the cube. The variation in δk0 seems to originate in the existence of a secondary component which dominates the outer parts.

In the text

Fig. C.1. Kinematic maps of M3G galaxies. From left to right, the mean velocity, velocity dispersion, h3 and h4 moments are shown. The values (in kilometres per second for the mean velocity and velocity dispersion) in lower-left corners correspond to maximum and minimum of the colour bars on the right of each plot. The red ellipse has the major axis length equal to Re and the median flattening of the galaxy within the effective radius (see Table 1). It shows the elliptical aperture used to extract global effective properties (e.g. σe, λRe, average values of kinemetry parameters, etc).

In the text

	Fig. C.1. Continued.
In the text

	Fig. C.1. Continued.
In the text

	Fig. C.1. Continued.
In the text

	Fig. C.1. Continued.
In the text