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

# A Deep View into the Nucleus of the Sagittarius Dwarf Spheroidal Galaxy with MUSE. II. Kinematic Characterization of the Stellar Populations

M. Alfaro-Cuello¹, N. Kacharov¹, N. Neumayer¹, P. Bianchini², A. Mastrobuono-Battisti¹,³, N. Lützgendorf⁴,

A. C. Seth⁵, T. Böker⁴, S. Kamann⁶, R. Leaman¹, L. L. Watkins⁷,⁸, and G. van de Ven⁸

Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

²Observatoire astronomique de Strasbourg, CNRS, UMR 7550, F-67000 Strasbourg, France

³Department of Astronomy and Theoretical Physics, Lund Observatory, Box 43, SE-221 00, Lund, Sweden

⁴European Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore, MD 21218, USA

⁵Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA

⁶Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK

⁷European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany

⁸Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17, A-1180 Wien, Austria

Received 2019 December 7; revised 2020 February 14; accepted 2020 February 17; published 2020 March 23

## Abstract

The Sagittarius dwarf spheroidal galaxy is in an advanced stage of disruption but still hosts its nuclear star cluster (NSC), M54, at its center. In this paper, we present a detailed kinematic characterization of the three stellar populations present in M54: young metal-rich (YMR); intermediate-age metal-rich (IMR); and old metal-poor (OMP), based on the spectra of ~6500 individual M54 member stars extracted from a large Multi-Unit Spectroscopic Explorer (MUSE)/Very Large Telescope data set. We find that the OMP population is slightly flattened with a low amount of rotation (~0.8 km s⁻¹) and with a velocity dispersion that follows a Plummer profile. The YMR population displays a high amount of rotation (~5 km s⁻¹) and a high degree of flattening, with a lower and flat velocity dispersion profile. The IMR population shows a high but flat velocity dispersion profile, with some degree of rotation (~2 km s⁻¹). We complement our MUSE data with information from Gaia DR2 and confirm that the stars from the OMP and YMR populations are comoving in 3D space, suggesting that they are dynamically bound. While dynamical evolutionary effects (e.g., energy equipartition) are able to explain the differences in velocity dispersion between the stellar populations, the strong differences in rotation indicate different formation paths for the populations, as supported by an N-body simulation tailored to emulate the YMR-OMP system. This study provides additional evidence for the M54 formation scenario proposed in our previous work, where this NSC formed via GC accretion (OMP) and in situ formation from gas accretion in a rotationally supported disk (YMR).

*Unified Astronomy Thesaurus concepts: Dwarf spheroidal galaxies (420); Galaxy nuclei (609); Globular star clusters (656); Sagittarius dwarf spheroidal galaxy (1423)*

## 1. Introduction

Nuclear star clusters (NSCs) are the densest stellar systems known in the universe with mass densities of ~$10^{6-7} M_\odot$ pc⁻³ (Walcher et al. 2005; Misgeld & Hilker 2011; Norris et al. 2014). They contain stellar masses from $10^5$ to $10^8 M_\odot$ in half-light radii of about 1–10 pc (Georgiev & Böker 2014), and are known to have extended and complex star-formation histories (e.g., Walcher et al. 2005; Kacharov et al. 2018). Nucleation is a common characteristic in galaxies across the Hubble sequence (Phillips et al. 1996; Carollo et al. 1998; Böker et al. 2002, 2004; Côté et al. 2006; Turner et al. 2012; Georgiev & Böker 2014), with similar occupation fractions for different galaxy cluster environments (Côté et al. 2006; Turner et al. 2012; den Brok et al. 2014; Muñoz et al. 2015; Sánchez-Janssen et al. 2019). The nucleation fraction peaks close to 100% in galaxies of $10^9 M_\odot$ and decreases toward both lower and higher galaxy masses (see Neumayer et al. 2020 for a recent review).

A large fraction of dwarf galaxies host low-mass NSCs that display characteristics similar to high-mass, metal complex globular clusters (GCs) in the Milky Way. This suggests that such GCs are former NSCs, remnants of the accretion of the host galaxy onto the Milky Way (Zinnecker et al. 1988; Böker 2008; Da Costa 2016). Cosmological merger simulations

show that Milky Way-like galaxies should harbor a non-negligible number of these GCs (Pfeffer et al. 2014; Kruijssen et al. 2019). This tidal stripping scenario has been strongly supported by the detection of supermassive black holes in more massive ultra-compact dwarf galaxies (e.g., Seth et al. 2014; Voggel et al. 2019).

Where detailed measurements of kinematics are possible, rotation seems to be a common dynamical ingredient of NSCs (Feldmeier et al. 2014; Nguyen et al. 2018) and is also commonly found in a high fraction of GCs (van de Ven et al. 2006; Bellazzini et al. 2012; Bianchini et al. 2013, 2018; Fabricius et al. 2014; Kacharov et al. 2014; Kimmig et al. 2015; Bellini et al. 2017; Kamann et al. 2018; Sollima et al. 2019). The presence of internal rotation can be an indicator of their formation mechanism (Mastrobuono-Battisti & Perets 2013, 2016; Hénault-Brunet et al. 2015; Gavagnin et al. 2016; Khoperskov et al. 2018; Mastrobuono-Battisti et al. 2019a) and can give important clues on their long-term dynamical evolution (e.g., Einsel & Spurzem 1999; Tiongco et al. 2018). Fabricius et al. (2014) and Kamann et al. (2018) found that GCs with high internal central rotation are more flattened, showing that internal rotation plays an important role in forming the shape of GCs.

NSC studies focus on galaxies with masses at or above the peak of the nucleation fraction ($10^9 M_\odot$), but no detailed works
are available for the low-mass galaxy regime. Thus, low-mass galaxy nucleation is not yet fully understood. In Alfaro-Cuello et al. (2019, hereafter Paper I), we started a study to understand NSCs in these unexplored low-mass galaxies and the chemodynamical properties they display. We looked for answers to these open questions by studying M54 (NGC 6715), the NSC of the Sagittarius dwarf spheroidal galaxy (Ibata et al. 1994)—hereafter Sgr dSph. This galaxy is currently being disrupted by the tidal field of the Milky Way (Ibata et al. 1997). The progenitor of the Sgr dSph galaxy is estimated to have a luminosity of ~$10^8 L_\odot$ ($M_v$ ~ -15; Niederste-Ostholt et al. 2012) and a dark halo mass with a lower limit of $6 \times 10^{10} M_\odot$ before infall (Gibbons et al. 2017; Mucciarelli et al. 2017). At a distance of 27.6 kpc (Sollima et al. 2009), M54 is located in the densest region of the Sagittarius stream at the photometric center of its host (Ibata et al. 1994; Mucciarelli et al. 2017).

In the literature, the central region of the Sgr dSph galaxy is described as two distinct populations: a metal-poor one with age of ~13 Gyr, identified as M54, and a metal-rich one with age of ~2 Gyr, identified as the nucleus of the Sgr dSph galaxy (Monaco et al. 2005; Siegel et al. 2007; Bellazzini et al. 2008). In Paper I, we suggested use of the term “M54” to describe the entire NSC of the Sgr dSph, as in the time of its discovery by Messier in 1778. We will follow the same approach in this paper, i.e., we will refer to M54 as the entire system and not just the metal-poor population.

In Paper I, we presented our large Multi-Unit Spectroscopic Explorer (MUSE) data set covering an area out to 2.5 times the half-light radius of M54 ($R_{HL}$ = 0.82, Harris 1996, 2010 edition), from which we extracted spectra of more than ~6600 member stars. From this extensive data set, we characterized the age and metallicity distributions of the stellar populations. We recovered the star-formation history of this NSC, detecting (at least) three stellar populations: young metal-rich (YMR), with a mean age of 2.2 Gyr and [Fe/H] = -0.04; intermediate-age metal-rich (IMR), with mean age of 4.3 Gyr and [Fe/H] = -0.29; and old metal-poor (OMP), with mean age of 12.2 Gyr and [Fe/H] = -1.41. From our findings, we suggested that the OMP population could be a merger remnant of two or more GCs driven to the center of the host galaxy by dynamical friction, thus explaining the large spread in both age and metallicity. We note that the OMP population sample might contain old, metal-poor stars that resided in the center of the Sgr dSph galaxy before the infall of the GCs that formed this massive dominant population we observe today. Thus, the large spread in age and metallicity that we measure might be partially due to metal-poor stars from the Sgr dSph. It is not possible to disentangle this with our data. The first encounter between the Sgr dSph galaxy and the Milky Way disk could have triggered the in situ star-formation episode (Tepper-García & Bland-Hawthorn 2018) of the YMR population in the center of the OMP population, resulting in a more flattened and centrally concentrated structure. The IMR population seems to be part of the Sgr dSph galaxy central field. However, the three different populations are interacting gravitationally, and thus, the IMR cannot be discarded from a dynamical analysis of the system.

Bellazzini et al. (2008) presented a kinematic analysis for a total of ~400 stars, including both metal-poor (OMP) and metal-rich (YMR + IMR) stars. The stars of both metallicity regimes were found to display consistent line-of-sight (LOS) velocities, as suggested by previous studies (e.g., Da Costa &

Armandroff 1995; Ibata et al. 1997; Monaco et al. 2005). According to these authors, the metal-rich stars follow a flat velocity dispersion profile of $\sigma$ = 10 km s⁻¹ within the central 9' radius. The metal-poor population follows a King profile (King 1966), which provides a good fit to the surface brightness, at least for the innermost region, as is usually found in GCs. This latter population has a maximum dispersion of $\sigma$ = 14.2 km s⁻¹ in the center and $\sigma$ = 5.3 km s⁻¹ at ~3.5 from the cluster center. Hence, the authors found differences in the kinematics between the two populations in a radius range of 1.5 < r < 6.5, suggesting different origins. In this picture, the metal-poor population enabled the collection of the gas that formed the metal-rich population, resulting in different velocity dispersion profiles. In addition, Bellazzini et al. (2008) found a weak signal of rotation for the metal-poor population < 2 km s⁻¹. This signal seems to increase at radii larger than 10', where a clear sign of rotation is displayed. Based on LOS velocities, Baumgardt & Hilker (2018) reported a central velocity dispersion for M54 of $\sigma_0$ = 16.2 km s⁻¹ (based on data from Bellazzini et al. 2008; Ibata et al. 2009; Carretta et al. 2010a; Watkins et al. 2015). Zooming into the center, kinematic observations and dynamical modeling from Ibata et al. (2009) and Baumgardt (2017) suggested the existence of an intermediate-mass black hole (IMBH) of $M_{BH}$ ~ $10^4 M_\odot$.

Using the exquisite observational data provided by MUSE, we are not just able to use integrated kinematics map of the NSC, as is usually done for extragalactic system. For M54, in contrast, we can measure to motion of individual stars, and hence, we can obtain the kinematics for each of the different stellar populations in M54. This provides valuable insights into the resulting complex kinematics of NSCs with multiple stellar structures. Since our MUSE data set covers 2.5 times the half-light radius of M54, we are able to evaluate the radial trends of the kinematic quantities up to significant distances from the cluster center.

In this second paper, we present a kinematic analysis of M54 to complement our previous study. We take into account the populations defined in Paper I, based on age and metallicity information. We present a brief description of the data in Section 2. We describe the method for the kinematic extraction in Section 3, including LOS velocity, velocity dispersion, and rotation. In Section 4, we present the kinematic results for M54 and each of its populations. We present a complement to our analysis using data from Gaia DR2 to reconstruct the 3D structure in Section 5. In Section 6, we show our development on ad hoc N-body simulations to support our interpretation of the data and of the formation history of M54. We discuss our findings in Section 7, conclude in Section 8, and briefly describe the future prospects in Section 9.

## 2. MUSE Data

As detailed in Paper I, the data set was acquired with the MUSE (Bacon et al. 2014), located at the UT4 of the Very Large Telescope at the Paranal Observatory in Chile as part of the programme 095.B-0585(A) (PI: Lützgendorf).

The 16 observed pointings (59.9 × 60.0 each) constitute a 4 × 4 mosaic covering a contiguous area out to ~2.5 times the half-light radius ($R_{HL}$ = 0.82, Harris 1996, 2010 edition) of M54. This is a total coverage of ~3.5 × 3.5, where 1' ≃ 8 pc at an assumed distance of 27.6 kpc (Sollima et al. 2009), or 28 × 28 pc. The MUSE data has a wavelength coverage of
<figure><img src="image_2.png" /><figcaption>Figure 1. Relation between the signal-to-noise logarithm and the LOS velocity errors obtained with UlySS.</figcaption></figure>

4800–9300 Å, with a mean spectral resolution of R ~ 2750 and spatial sampling of 0.2'2 pix⁻¹.

For more details about the data, single stellar extraction, and membership estimates, see Paper I.

### 3. Analysis

In this section, we describe the methods we use for the kinematic extraction that are applied to the samples presented in Section 4.

#### 3.1. LOS Velocity

In Paper I, we described the determination of physical parameters and LOS velocities of the stars using University of Lyon Spectroscopic Analysis Software (ULySS; Koleva et al. 2009). To obtain the stellar atmospheric parameters, ULySS interpolates and fits to the observed spectra a stellar library of synthetic spectral templates, characterized by metallicity [Fe/H], surface gravity, and temperature. We use the synthetic spectroscopic grid, built on the basis of the ELODIE 3.2 library (Wu et al. 2011). We present in Figure 1 the relation between the signal-to-noise logarithm of the stars and the error in the velocity measurements obtained with UlySS. We exclude from the sample all stars with error > 20 km s⁻¹ leaving 7123 stars. This number includes all extracted stars, and we describe the membership cleaning of the sample in Section 3.2. The median of the uncertainties of the LOS velocities estimated by UlySS is 2.18 km s⁻¹. To test the robustness of our results, we perform the same kinematic analysis for stars with errors lower than 5 km s⁻¹. This subsample includes 5446 stars and we obtain consistent results. We corrected the LOS velocity of the stars for the effect of perspective rotation using Equation 6 in van de Ven et al. (2006).

#### 3.2. Rotation and Velocity Dispersion

Constraining the rotation of the populations of M54 can provide fundamental constraints on their origin. We measure the rotation and velocity dispersion of M54 simultaneously, using a discrete Bayesian approach similar to the procedure described in Cordero et al. (2017) and Koch et al. (2018) where we consider the LOS velocities of the stars and their respective errors. We approximate our velocity dispersion profile with a

Plummer model (Plummer 1911):

$$ \sigma(r)^2 = \frac{\sigma_0^2}{\sqrt{1 + \frac{r^2}{a^2}}}, \qquad (1) $$

where σ₀ and a are the central velocity dispersion and the Plummer radius, respectively, both treated as free parameters. The mean LOS velocity (<V<sub>LOS</sub>) is also a free parameter. We approximate the rotation profile with an analytical model (Mackey et al. 2013; Kacharov et al. 2014; Cordero et al. 2017):

$$ V_{\text{rot}} \sin i = \frac{2V_{\text{max}}}{r_{\text{peak}}} \times \frac{X_{\text{PA}}}{1 + (X_{\text{PA}}/r_{\text{peak}})^2}, \qquad (2) $$

where V<sub>rot</sub> sin i is the rotation amplitude at position X<sub>PA</sub> within a factor of the cluster’s inclination (i) with respect to the LOS. The independent variable X<sub>PA</sub> represents the distance of the different points from the center of M54 along the equatorial axis. For this model, we set three free parameters: V<sub>max</sub>, the maximum rotation amplitude at the projected radius r<sub>peak</sub>, and the rotation axis position angle (PA; measured from north 0° to east 90°). We set a weak Gaussian prior on r<sub>peak</sub> centered at 0 with a standard deviation of 4’—about five times the half-light radius and ~0.5 times the tidal radius of the cluster. We use a Gaussian mixture likelihood model to exclude foreground stars, where the mixing factor is the membership fraction, from which we calculate the membership probability of each star. To summarize, we account for a total of seven free parameters in our discrete kinematic model: σ₀, a, <V<sub>LOS</sub>>, V<sub>max</sub>, r<sub>peak</sub>, PA, and membership fraction. We optimize the Gaussian likelihood function (see Equation (2) and (5) in Cordero et al. 2017) using the affine-invariant Markov chain Monte Carlo (MCMC) algorithm EMCEE (Goodman & Weare 2010; Foreman-Mackey et al. 2013). We use 500 walkers and 5000 iterations to assure all parameters converge. We use this method since evaluating discrete data provides more accurate results than other methods that require binning the data, which can end in loss of information.

To be consistent with other studies, we measure the average rotation amplitude (A<sub>rot</sub>) as the half of the difference between the maximum likelihood mean velocity on the two sides of the rotation axis. For the analysis in Section 4.3, we additionally measure the rotation at the half-light radius (A<sub>HL</sub>) and its uncertainty, as the median and standard deviation of the rotation model posterior distribution at that radius.

### 4. Kinematic Extraction

In this section, we present the results of the kinematic extraction applied to: (i) all M54’s member stars as a whole— with no distinction in populations—and (ii) the three populations in M54 identified in Paper I.

The results illustrated in this section are summarized in Table 1. For an easier comparison, we include in this table some of the population parameters from Paper I (i.e., [Fe/H], age, ellipticity, etc.).
<table><caption>**Table 1**
Summary of the Observed Properties of M54 and Its Stellar Populations</caption><thead><tr><th>Populations</th><th>M54 members</th><th>YMR</th><th>IMR</th><th>OMP</th></tr></thead><tbody><tr><td>[Fe/H]</td><td>-2.5 to 0.5</td><td>-0.04 ± 0.01<sup>a</sup></td><td>-0.29 ± 0.01<sup>a</sup></td><td>-1.41 ± 0.01<sup>a</sup></td></tr><tr><td>σ<sub>[Fe/H]</sub></td><td>...</td><td>0.12 ± 0.01<sup>a</sup></td><td>0.16 ± 0.01<sup>a</sup></td><td>0.24 ± 0.01<sup>a</sup></td></tr><tr><td>Age (Gyr)<sup>a</sup></td><td>0.5 to 14</td><td>2.16 ± 0.03<sup>a</sup></td><td>4.28 ± 0.09<sup>a</sup></td><td>12.16 ± 0.05<sup>a</sup></td></tr><tr><td>σ<sub>Age</sub> (Gyr)</td><td>...</td><td>0.20 ± 0.03<sup>a</sup></td><td>1.16 ± 0.07<sup>a</sup></td><td>0.92 ± 0.04<sup>a</sup></td></tr><tr><td>ε</td><td>0.13 ± 0.03</td><td>0.31 ± 0.10<sup>a</sup></td><td><sup>b</sup></td><td>0.16 ± 0.06<sup>a</sup></td></tr><tr><td>ε<sub>HL</sub></td><td>0.08 ± 0.04</td><td>0.15 ± 0.13</td><td><sup>b</sup></td><td>0.09 ± 0.06</td></tr><tr><td>PA (deg)</td><td>8.83 ± 6.62</td><td>4.56 ± 10.84<sup>a</sup></td><td><sup>b</sup></td><td>16.83 ± 12.81<sup>a</sup></td></tr><tr><td>Median V<sub>LOS</sub> (km s<sup>-1</sup>)</td><td>141.54 ± 0.21</td><td>142.01 ± 0.52</td><td>142.61 ± 0.59</td><td>141.22 ± 0.27</td></tr><tr><td>σ<sub>0</sub> (km s<sup>-1</sup>)</td><td>15.91 ± 0.33</td><td>12.56 ± 0.72</td><td>15.21 ± 0.89</td><td>15.30 ± 0.54</td></tr><tr><td>A<sub>rot</sub> (km s<sup>-1</sup>)</td><td>1.16 ± 0.36</td><td>4.91 ± 0.90</td><td>2.26 ± 1.24</td><td>0.75 ± 0.55</td></tr><tr><td>A<sub>HL</sub> (km s<sup>-1</sup>)</td><td>1.60 ± 0.29</td><td>7.30 ± 0.84</td><td>2.00 ± 0.89</td><td>0.73 ± 0.38</td></tr><tr><td>V<sub>max</sub> (km s<sup>-1</sup>)</td><td>1.88 ± 0.70</td><td>7.98 ± 1.01</td><td>3.15 ± 1.95</td><td>1.58 ± 1.14</td></tr><tr><td>Rotation axis (deg)</td><td>2.9 ± 8.8</td><td>-0.1 ± 7.2</td><td>-5.9 ± 26.4</td><td>4.3 ± 23.7</td></tr><tr><td>Number of stars</td><td>6537</td><td>630</td><td>649</td><td>2550</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Values from Paper I.

<sup>b</sup> Measurements that did not converge to a value. σ<sub>[Fe/H]</sub> and σ<sub>Age</sub> correspond to the [Fe/H] and age intrinsic spreads, respectively. ε: ellipticity estimated at the entire field of view (FOV). ε<sub>HL</sub>: ellipticity estimated at the half-light radius. V<sub>LOS</sub>: LOS velocity. σ<sub>0</sub>: central velocity dispersion. A<sub>rot</sub>: average rotation amplitude. A<sub>HL</sub>: rotation at the half-light radius. V<sub>max</sub>: maximum rotation. Note: The number of stars in the M54 sample is larger than the sum of the three populations because it includes horizontal branch stars (which were excluded for the age estimates) and stars that could not clearly be attributed to either of the three subpopulations when we created the multi-Gaussian model in Paper I.

## 4.1. Kinematics of All M54 Member Stars

After excluding background stars from our stellar sample as described in Section 3.2, the final sample of M54 member stars includes 6537 stars. These stars occupy a large range in age (0.5–14 Gyr) and metallicity (−2.5 < [Fe/H] < 0.5). The use of this overall sample mimics the kinematic extraction from a single stellar system, as is commonly the case for more distant NSCs, for which it is impossible to resolve individual stars or to separate different stellar populations.

### 4.1.1. LOS Velocity and Velocity Dispersion

We include the LOS velocity map for the resolved stars in the top left panel of Figure 2. The dashed black circle shows the half-light radius of M54 ($R_{HL}$ = 0.82, Harris 1996, 2010 edition). For illustration purposes, in the top panel of Figure 2, we present Voronoi-binned maps for the mean LOS velocity (middle) and velocity dispersion (right) for our observed M54 member stars. We use the Voronoi tessellation code by Cappellari & Copin (2003) as described by Kamann et al. (2018). We aim for ~45 stars per bin.

Based on our discrete model, the mean V<sub>LOS</sub> value is 141.54 ± 0.21 km s<sup>-1</sup>. We include in the bottom right panel the best-fitting Plummer profile to the data (solid cyan line), and the ±3σ uncertainty (dashed cyan lines). The cyan circles show the velocity dispersion from the data points estimated in radial bins of 0.3, and the vertical dashed line shows the half-light radius of M54. The median of the best-fit central velocity dispersion is σ<sub>0</sub> = 15.91 ± 0.33 km s<sup>-1</sup>.

### 4.1.2. Rotation

The results of the rotation analysis for M54 are presented in the bottom left and middle panels of Figure 2. The left panel shows the velocity gradient with respect to a line perpendicular to the rotation axis. The solid black line is the median of the best-fit rotation models. The dashed black lines show ±3 times the velocity dispersion. In the bottom middle panel of Figure 2,

we show the best-fit rotation models to the discrete data obtained using Equation (2), as a solid cyan line, and the ±3σ uncertainty as dashed cyan lines. For representation, the cyan circles represent the rotation profile derived as the difference between the median velocity and the systemic velocity for overlapping bins of 0.6 along the line perpendicular to the rotation axis. The horizontal and vertical error bars represent the radial bin size and the uncertainties in the offset of the median velocity, respectively. The vertical dashed lines show the half-light radius of M54 ($R_{HL}$ = 0.82, Harris 1996, 2010 edition). The best-fit rotation axis corresponds to 2.9 ± 8.8. We detect a low—but still considerable—amount of rotation in M54. We obtain a rotation amplitude of A<sub>rot</sub> = 1.16 ± 0.36 km s<sup>-1</sup>, and a maximum rotation of V<sub>max</sub> = 1.88 ± 0.70 km s<sup>-1</sup>.

## 4.2. Kinematics of M54's Populations

In this section, we study the kinematics of each of the three populations identified via age and metallicity in Paper I: YMR (630 stars), IMR (649 stars), and OMP (2550 stars). For better kinematic characterization, we do not exclude stars with ages with relative errors >40%, thus, the number of stars for each population is slightly larger than those considered in Paper I. The number of stars in the full M54 sample is larger than the sum of stars in the three populations, because it also includes horizontal branch stars (~650) and stars that could not clearly be attributed to either of the three subpopulations when we created the multi-Gaussian model in Paper I.

### 4.2.1. LOS Velocity and Velocity Dispersion

In the left, middle, and right panels of Figure 3, we present the LOS velocity maps of the resolved stars, and for representation, the Voronoi-binned velocity and velocity dispersion maps for each of the populations in M54: YMR, IMR, and OMP, from top to bottom, respectively. For the Voronoi maps, we set a target number of stars per bin
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Kinematics for all the 6537 member stars of M54. Top panels: LOS velocity map from resolved stars (left), Voronoi-binned maps for the mean LOS velocity (middle) and velocity dispersion (right). The overplotted solid gray lines show the median of the best-fit rotation axis. The overplotted dashed black circle in the left panel shows the half-light radius of M54 ($R_{HL} = 0.82$, Harris 1996). The green ellipse in the top middle and right panels is overplotted using the half-light radius, PA and ellipticity estimated in Paper I fitting a King profile. Bottom left panel: velocity gradient with respect to a line perpendicular to the rotation axis. The solid black line shows the median of the best-fit rotation models. The dashed black lines show ±3 times the velocity dispersion. Bottom middle panel: rotation profile for M54 member stars (solid cyan line) obtained as the median of the best-fit rotation models using Equation (2). The cyan circles represent the rotation profile derived as the difference between the median velocity and the systemic velocity for overlapping bins of 0.6 along the line perpendicular to the rotation axis. The horizontal and vertical error bars represent the radial bin size and the uncertainties in the offset of the median velocity, respectively. Bottom right panel: velocity dispersion profile for the M54 members obtained as the median of the best-fit Plummer profiles (solid cyan line). The dashed cyan lines show the ±3σ uncertainty for each best-fit model. The cyan circles show the velocity dispersion from the data points estimated in radial bins of 0.3. The vertical dashed lines in both panels show the half-light radius of M54 ($R_{HL} = 0.82$, Harris 1996, 2010 edition).</figcaption></figure>

depending of the size of the sample of the populations: 25 for the OMP, and seven for the IMR and YMR populations.

We found very similar mean LOS velocities for the three populations. The mean velocities are 142.01 ± 0.52 km s⁻¹ for the YMR, 142.61 ± 0.59 km s⁻¹ for the IMR, and 141.22 ± 0.27 km s⁻¹ for the OMP. The LOS velocities are consistent when only considering stars with errors lower than 5 km s⁻¹, showing variations of ~0.3 km s⁻¹. From the Voronoi-binned LOS velocity maps, we observe a clear velocity gradient for the YMR population.

We present, in the right panels of Figure 4, the velocity dispersion profile for each population: YMR (top), IMR (middle), and OMP (bottom). The profile is obtained as the median of the best-fit Plummer models to the discrete data. The dashed lines show the ±3σ uncertainty for each profile. For representation, in each case, the colored circles show the velocity dispersion from the data points estimated in radial bins of 0.3. The vertical dashed lines in both panels show the half-light radius of M54 ($R_{HL} = 0.82$, Harris 1996, 2010 edition). We estimate a central velocity dispersion of

$\sigma_0 = 12.56 \pm 0.72 \text{ km s}^{-1}$ for the YMR, and $\sigma_0 = 15.21 \pm 0.89 \text{ km s}^{-1}$ for the IMR. For the YMR and IMR velocity dispersion profiles, we observe that they rather follow a flat distribution along the observed radius; thus, the best-fit Plummer profiles do not provide a good description of the data for these two populations. Consistently, the best-fit velocity dispersion profiles for the YMR and IMR populations are flat when we only consider stars with errors lower than 5 km s⁻¹, displaying slightly flatter profiles with central σ values lower by ~1 km s⁻¹. In contrast, the OMP velocity dispersion profile (blue) is very well fit by a Plummer model. The central velocity dispersion for the OMP population is $\sigma_0 = 15.30 \pm 0.54 \text{ km s}^{-1}$.

### 4.2.2. Rotation

In the left panels of Figure 4, we present the velocity gradient with respect to a line perpendicular to the rotation axis for the three populations: YMR (top), IMR (middle), and OMP (bottom). The solid black line shows the median of the best-fit rotation models. The dashed black lines show ±3 times the
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Kinematic extraction for the three M54 populations: YMR in top panels, IMR in middle panels, and OMP in bottom panels. The left panel shows the LOS velocity map for the stars in each population. The middle and right panels show the Voronoi-binned maps for the mean LOS velocity and velocity dispersion, respectively. Overplotted dashed black circles in the left panels show the half-light radius of M54 ($R_{HL} = 0.82$, Harris 1996). The green ellipses in the middle and right panels are overplotted using the half-light radius, PA and ellipticity estimated in Paper I for a King profile. Since these values were not obtained for the IMR population, we overplotted a spherical shape. The overplotted gray solid lines show the median of the best-fit rotation axis.</figcaption></figure>

velocity dispersion. In the middle panels, we present the median of the best-fit rotation profiles obtained using Equation (2) as solid lines for the three populations. The dashed lines of the same colors represent the ±3σ uncertainties. The vertical dashed lines show the half-light radius of M54 ($R_{HL} = 0.82$, Harris 1996, 2010 edition). As an illustration, the colored circles represent the rotation profile derived as the difference between the median velocity and the systemic velocity for overlapping bins of 0.6 along the line perpendicular to the rotation axis. The horizontal and vertical error bars represent the radial bin size and the uncertainties in the offset of the median velocity, respectively.

The highest maximum rotation is displayed by the YMR population with $V_{max} = 7.98 \pm 1.01$ km s⁻¹, with a median rotation axis at −0.1 ± 7.2. The second highest maximum rotation is displayed by the IMR population with a median of $V_{max} = 3.15 \pm 1.95$ km s⁻¹, less than half that of the YMR population, and a median rotation axis at −5.9 ± 26.4. The OMP population shows a low amount of rotation, a best-fit value of maximum rotation of $V_{max} = 1.58 \pm 1.14$ km s⁻¹, with a rotation axis at 4.3 ± 23.7. In order to compare with other studies, we also compute the average amplitudes of
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left panels: velocity gradient with respect to a line perpendicular to the rotation axis for the three populations: YMR (top), IMR (middle), and OMP (bottom). The solid black line shows the median of the best-fit rotation models and the shaded area corresponds to ±3σ. The dashed black lines show ±3 times the velocity dispersion. Middle panels: rotation profile for the three populations obtained as the median of the best-fit rotation models using Equation (2). The dashed lines show the ±3σ. The colored circles represent the rotation profile derived as the difference between the median velocity and the systemic velocity for overlapping bins of 0.6′ along the line perpendicular to the rotation axis. The horizontal and vertical error bars represent the radial bin size and the uncertainties in the offset of the median velocity, respectively. Right panels: velocity dispersion profile for the three populations obtained as the median of the best-fit Plummer profiles (solid lines). In all panels, the vertical dashed lines show the half-light radius of M54 ($R_{HL}$ = 0.82, Harris 1996, 2010 edition).</figcaption></figure>

rotation<fnref n="9" /> of the three populations as follows: $A_{\text{rot(OMP)} = 0.75 \pm 0.55 \text{ km s}^{-1}$, $A_{\text{rot(IMR)} = 2.26 \pm 1.24 \text{ km s}^{-1}$, and the fastest rotator with $A_{\text{rot(YMR)} = 4.91 \pm 0.90 \text{ km s}^{-1}$. The differences in the velocity dispersion and rotation properties suggest different origins for these populations.

### 4.3. $V/\sigma_0$ versus $\epsilon$

The ratio between the rotation and central velocity dispersion versus ellipticity ($V/\sigma_0, \epsilon$) diagrams have been introduced and are widely used to evaluate how the rotation affects the shape of galaxies (e.g., Binney 2005; Cappellari et al. 2007; Emsellem et al. 2011). However, in the last few years, these diagrams have also been applied to GCs (e.g., Bellazzini et al. 2012; Bianchini et al. 2013, 2018; Fabricius et al. 2014; Kacharov et al. 2014; Kimmig et al. 2015; Lardo et al. 2015; Kamann et al. 2018). In this approach, several difficulties arise for both clusters and galaxies, because: (i) the rotation

measured in clusters depends on the inclination angle of the rotation axis, which is measurable only when 3D kinematics are available (LOS velocities and proper motions, e.g., Bianchini et al. 2018; Sollima et al. 2019); (ii) both ellipticity and rotation change as a function of radius (e.g., Geyer et al. 1983); and (iii) the anisotropy also varies with radius (e.g., van de Ven et al. 2006; Jindal et al. 2019). In spite of these issues, ($V/\sigma_0, \epsilon$) diagrams still provide a good first diagnostic of the role of rotation on the shape of these stellar systems, e.g., how flattened they are (Fabricius et al. 2014; Kamann et al. 2018). In Figure 5, we present a ($V/\sigma_0, \epsilon$) diagram where we add the measurements for the whole M54 sample (gray), the YMR (red), and OMP (blue) populations. We do not include the IMR population in this section, since the method used to obtain the 2D morphology parameters did not converge to an ellipticity value. Since different definitions of $V/\sigma_0$ are used in the literature, we decide to compute the $V/\sigma_0$ in two different ways: (i) based on the rotation and ellipticity at the half-light radius ($V_{\text{HL}}/\sigma_0, \epsilon_{\text{HL}}$) and represented by circles, and (ii) based on the rotation amplitude and ellipticity over the entire FOV ($A_{\text{rot}}/\sigma_0, \epsilon$), represented by squares. The ellipticity at the half-

<footnote n="9">⁹ As half the difference of the mean velocities on the two sides of the rotation axis.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** ($V/\sigma_0$, $\epsilon$) diagram. We present our estimates for the entire sample of stars (M54), the YMR and OMP populations in gray, red, and blue, respectively. The circles show the estimates based on the rotation and ellipticity at the half-light radius ($V_{HL}/\sigma_0$, $\epsilon_{HL}$), the squares, on the rotation amplitude and ellipticity over the entire FoV ($A_{rot}/\sigma_0$, $\epsilon$). For comparison, we include values for GCs from the literature (e.g., Bellazzini et al. 2012; Bianchini et al. 2013; Kacharov et al. 2014) and in green squares, the GC sample from Kamann et al. (2018). The solid black line gives the relation for an isotropic oblate rotator as Binney (2005). The dashed black lines show the relation for oblate rotators with anisotropy parameters of $\delta = 0.05, 0.1, 0.15, 0.2$, from left to right, as Cappellari et al. (2007).</figcaption></figure>

light radius is derived using the method described in Paper I used to estimate the ellipticity across the entire FoV. We explore the possibility that the inclination angle, unknown a priori, could be different from an edge-on view and, thus, introduce a correction for inclination as described in Cappellari et al. (2007). We present the respective corrected values in Table 2. We will show in Section 5 that by combining LOS data and Gaia DR2 proper motions, we can estimate that the most likely inclination angle is ~90° (edge-on view). However, inclination angles within 60° and 90° cannot be excluded within the $1\sigma$ uncertainties. Therefore, we illustrate the effects of a 60° inclination angle, and consider the corrected value as a conservative (upper limit) estimate for the intrinsic ($V/\sigma_0$, $\epsilon$).

In Figure 5, we include for comparison measurements for GCs from Bellazzini et al. (2012), Bianchini et al. (2013), Kacharov et al. (2014), Kamann et al. (2018), who consistently measure ($V/\sigma_0$) at the half-light radius or beyond. We overplot as a solid black line the relation for an isotropic oblate rotator as in Binney (2005), and the dashed black lines show the relation for oblate rotators with the global anisotropy $\delta$ (see Binney & Tremaine 1987) parameters $\delta = 0.05, 0.1, 0.15, 0.2$ from left to right, as presented in Cappellari et al. (2007).

For the OMP population, we observe a low amount of rotation, with only minor differences when estimating this value at the half-light radius $R_{HL}$ (blue circle) or over the entire FOV of the sample (blue square). In spite of a low amount of rotation, the ellipticity estimated over the entire FOV shows that this population is flattened, similar to fast rotating GCs such as the cases plotted in Figure 5: $\omega$ Cen (black diamond) and M22 (black pentagon). We will discuss a possible explanation for this in Section 7.3.

<table><caption>**Table 2**
Corrected ($V\sigma$) and Ellipticity ($\epsilon$) Values</caption><thead><tr><th></th><th>M54</th><th>YMR</th><th>OMP</th></tr></thead><tbody><tr><td>$A_{rot}/\sigma_0$</td><td>0.08</td><td>0.44</td><td>0.06</td></tr><tr><td>$\epsilon$</td><td>0.18</td><td>0.45</td><td>0.22</td></tr><tr><td>$V_{HL}/\sigma_0$</td><td>0.11</td><td>0.70</td><td>0.05</td></tr><tr><td>$\epsilon_{HL}$</td><td>0.11</td><td>0.21</td><td>0.12</td></tr></tbody></table>

**Note.** All values were obtained correcting for inclination as described in Cappellari et al. (2007).

Our rotation estimates for M54 and its OMP population are consistent with the estimate for the metal-poor population by Bellazzini et al. (2012), who considered an ellipticity of $\epsilon = 0.06$ from Harris (1996), 2010 edition. A more consistent agreement between the rotation and ellipticity of the OMP population reported by Bellazzini et al. (2012) was expected since their sample does not include young and metal-rich stars.

### 4.4. Energy Equipartition

In this section, we explore the possibility that differences in the kinematics could arise from dynamical relaxation processes, namely energy equipartition, and not necessarily only from different formation mechanisms of the different populations.

Dynamical interactions between stars of different masses naturally bring a stellar system into a mass-segregated configuration, with massive stars losing energy and sinking to the center. In parallel, stars with different masses will display different kinematics, a process known as partial energy equipartition (see e.g., Spitzer 1969; Trenti & van der Marel 2013; Bianchini et al. 2016) that acts on the timescale of a few relaxation times. In particular, massive stars are expected to have a lower velocity dispersion than low-mass stars, with a velocity dispersion ratio that depends on the relaxation state of the cluster and the stellar mass difference between the populations. In this study, the YMR population is characterized by an average stellar mass$^{10}$ of $1.53 M_{\odot}$, while that of the OMP population is a factor ~2 smaller, $0.82 M_{\odot}$. Therefore, different kinematics could arise during their ~2 Gyr coevolution (i.e., the age of the YMR population).

Following Equation (3) of Bianchini et al. (2016), the ratio between the velocity dispersions of two stellar components due to partial energy equipartition is given by:

$$ \frac{\sigma_{\text{OMP}}}{\sigma_{\text{YMR}}} = \exp\left(-\frac{1}{2}\frac{m_{\text{OMP}} - m_{\text{YMR}}}{m_{\text{eq}}}\right) \quad (3) $$

with $m_{YMR}$ and $m_{OMP}$ the typical masses of the YMR and OMP component and $m_{eq}$ the equipartition parameter quantifying what degree of energy equipartition the system was able to reach (i.e., how strong the kinematic differences are expected to be). While the equipartition parameter is not known a priori, it strongly depends on the relaxation state of a stellar system: clusters that have lived many relaxation times are characterized by a higher degree of energy equipartition.

We can use Equation (6) of Bianchini et al. (2016) to estimate the value of $m_{eq}$, starting from the current core

<footnote>$^{10}$ We obtained the stellar mass for each star from the best-fit isochrone as described in Paper I.</footnote>
relaxation time of M54 given by Harris (1996), $T_{rc} = 1.738 \times 10^8$ yr. This implies that M54 has dynamically evolved for ~12 relaxation times, which yields $m_{eq} = 2.03 M_{\odot}$. Using this value in Equation (3), we estimate that dynamical processes can produce a ratio between the velocity dispersion of the YMR and OMP population of $\sigma_{OMP}/\sigma_{YMR} = 1.19$. This difference is fully consistent with the one we observe (~1.22).

In summary, the effect of partial energy equipartition due to dynamical relaxation can indeed explain the kinematic differences between the YMR and OMP population. In other words, the two populations need not necessarily have formed with different velocity dispersions in order to display the currently observed differences. We caution, however, that further studies of the interplay between energy equipartition and rotation are needed for a comprehensive interpretation of the dynamical evolution of the kinematic differences.

## 5. M54 kinematics with Gaia DR2

The diagnostic diagram ($V/\sigma_0$, $\epsilon$) suffers from the degeneracy of information of the inclination angle, i.e., unknown intrinsic amount of angular momentum. This limitation can be compensated for by using 3D velocity data from Gaia DR2. The recent availability of precision astrometry from Gaia Data Release 2 (DR2; Gaia Collaboration et al. 2016, 2018) opens the possibility of studying the kinematics of GCs in 3D (e.g., Sollima et al. 2019). Here, we exploit for the first time the available Gaia proper motions with the goal of characterizing the intrinsic dynamical properties of M54's populations.

We cross match our MUSE M54 stellar sample with the Gaia DR2 catalog (Gaia Collaboration et al. 2016, 2018) finding a total of 638 OMP stars, 109 YMR stars, and 32 IMR stars, for which the full 3D velocity vector is measured. These stars have magnitudes $G < 20$, and are within ~2'0 from M54's center. However, M54 represents a challenging environment for Gaia astrometry due to its relatively large distance and the resulting high degree of crowding. To obtain reliable measurements, we need to further restrict the sample to only stars with high-quality proper motion measurements. For this purpose, we perform quality cuts, following some of the procedures illustrated in Lindegren et al. (2018) and Vasiliev (2019), to eliminate stars with bad astrometric measurements and those strongly affected by crowding. These cuts, based on parameters provided in the Gaia DR2 catalog,<fnref n="11" /> include:

1. astrometric_gof_al < 0.5

2. astrometric_excess_noise < 1

3. unit weight error, $uw = < 1.2$ (see Lindegren et al. 2018 for definition)

4. phot_bp_rp_excess_factor < 2.0 + 0.06(bp-rp)², with bp-rp the color in Gaia filters.

5. proper motions errors < 0.5 mas yr⁻¹ (corresponding to ~60 km s⁻¹ at M54's distance).

Our final sample consists of a total of 108 OMP stars and 15 from the YMR population, with an average uncertainty in proper motion measurements of 0.19 mas yr⁻¹ (corresponding to ~25 km s⁻¹). The final sample for the IMR populations consist of two stars only. Since this is not sufficient to perform

the analysis, we will just consider the YMR and OMP populations.

### 5.1. Rotation with Gaia DR2

We transform the positions and velocities from celestial to Cartesian coordinates using Equation (2) of Gaia Collaboration et al. (2018) and Equation (1) of van de Ven et al. (2006; see also Bianchini et al. 2018). We correct the LOS velocities and proper motions for perspective rotation following Equation (6) in van de Ven et al. (2006). Finally, when converting proper motions from mas yr⁻¹ to km s⁻¹, we assume a distance of 27.6 kpc (Sollima et al. 2009). As a further test, we repeat the analysis using a distance of 26.5 kpc (Harris 1996, 2010 edition) and 28.4 kpc (Siegel et al. 2011), obtaining consistent results.

We measure the mean motions of the two proper motion components using the likelihood employed in Bianchini et al. (2018, Equations (2) and (3)). For the OMP stars we obtain mean motions of ($\mu_x, \mu_y$) = (2.80 ± 0.03, 1.40 ± 0.02) mas yr⁻¹, and ($\mu_x, \mu_y$) = (2.77 ± 0.09, 1.48 ± 0.08) mas yr⁻¹ for the YMR stars. These values are consistent with each other and with the value reported by Vasiliev (2019), who made no distinction between populations in M54. Together with the fact that the mean LOS velocities of the OMP and YMR components are also consistent with each other, this further indicates that the two stellar populations are comoving, and therefore, they belong to the same bound stellar system.

To search for a signature of rotation, we consider polar coordinates in the plane of the sky and analyze the tangential component of proper motions ($\mu_t$). Both populations show a mean value of $\mu_t$ consistent with zero ($\mu_{t, OMP} = 0.012 \pm 0.030$ mas yr⁻¹ and $\mu_{t, YMR} = -0.058 \pm 0.113$ mas yr⁻¹), indicating no signal of rotation, within their $1\sigma$ uncertainties. However, we note that the sample size of the YMR population is composed of only 15 stars and, moreover, the putative rotation signal that we are trying to measure ($< 10$ km s⁻¹) is below the nominal systematic uncertainties for Gaia DR2 data (0.07 mas yr⁻¹, Lindegren et al. 2018). Therefore, a presence of rotation cannot be excluded with this analysis and the current data.

### 5.2. Rotation from 3D Kinematics

Since our sample consists of the full three-dimensional velocity vectors, we can estimate the intrinsic rotation exploiting simultaneously the three velocity components, following the likelihood method described by Sollima et al. (2019), likelihood in Equation (3). We assume a constant rotation amplitude<fnref n="12" /> within the cluster and take into consideration the discrete velocity measurements, their uncertainties and their covariance matrix. We sample the likelihood and derive the $1\sigma$ error using emcee by Foreman-Mackey et al. (2013) keeping as free parameters the PA of the rotation axis in the plane of the sky ($\theta_0$, measured from west of north), the inclination angle of the rotation axis with respect to the LOS ($i$, with $i = \pi/2$, corresponding to an edge-on view), and the amplitude of the intrinsic rotation ($A_{\text{rot}}$, measured in km s⁻¹). We fix a velocity dispersion of 13 km s⁻¹ (as derived globally from LOS velocities only) for the modeling. We repeated the

<footnote n="11">¹¹ For a detailed description of the Gaia parameters, see Lindegren et al. (2018).</footnote>

<footnote n="12">¹² This is a first-order approximation suitable for the low number of data points available.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Results of the MCMC sampling of intrinsic 3D rotation of the OMP component (left panel; 108 data points) and of the YMR component (right panel; 15 data points). The sampled parameters are the PA of the rotation axis in the plane of the sky ($\theta_0$), the amplitude of the rotation ($A_{\text{rot}}$), and the inclination angle of the rotation axis with respect to the LOS ($i$). No rotation is observed for the OMP stars, while clear rotation, mostly along the LOS velocity component, is measured for the YMR stars. The blue lines show the mean parameter value of the respective axis.</figcaption></figure>

analysis using different values between 10 and 15 km s⁻¹ and saw no differences in the final results.

The results obtained for the OMP and YMR components are shown in Figure 6 and reported in Table 3. We do not detect intrinsic rotation in the OMP component (2.0⁺¹.⁴⁻¹.⁸ km s⁻¹). On the other hand, the YMR population shows a strong signature of rotation (13.5⁺⁶.⁴⁻⁵.⁸ km s⁻¹), predominantly along the LOS (note that this signal is dominated by low-number statistics, as reflected by the large statistical uncertainties), since the recovered inclination angle $i$ is consistent with an edge-on view. In Figure 7, we show the result of the simultaneous modeling of the three velocity components for the YMR component.

## 6. Comparison with N-body Simulations.

To understand the origin of the morphological and kinematical structure of the system, we simulated a two-component M54-like cluster with a total mass equal to 1.4 × 10⁶ M☉. The modeled cluster is characterized by a YMR/OMP mass ratio of 0.20 (based on the number ratio of stars), as suggested by the observations. The OMP population is represented by a non-rotating, spherical King (1966) profile with a total mass of 1.17 × 10⁶ M☉, $W_0 = 8.6$, and a core radius of 0.72 pc (Harris 1996, 2010 edition). The remaining mass of the cluster is included as a flattened and centrally concentrated YMR population (50% of the mass is enclosed within 3 pc from the center of the cluster) that formed in situ from gas captured by the gravitational potential of the OMP component. The gas kept its angular momentum and formed a rotating stellar component. This disk-like component has an ellipticity ranging from $\epsilon \approx 0.1$ in the central region to $\epsilon = 0.5–0.6$ at radii larger than 6 pc. The YMR population initially rotates with a peak velocity of 12 km s⁻¹, a value similar to the maximum rotational velocity of the YMR population observed in M54.

<figure><img src="image_8.png" /><figcaption>**Table 3**

3D Modeling of the Intrinsic Rotation Using the Three Velocity Components</figcaption></figure>

<table><thead><tr><th>Population</th><th>θ<sub>0</sub> (deg)</th><th>A<sub>rot</sub> (km s<sup>-1</sup>)</th><th>i (deg)</th></tr></thead><tbody><tr><td>YMR</td><td>32.4<sup>+24.6</sup><sub>-18.5</sub></td><td>13.0<sup>+6.9</sup><sub>-6.6</sub></td><td>77.9<sup>+9.1</sup><sub>-18.9</sub></td></tr><tr><td>OMP</td><td>-17.0<sup>+96.9</sup><sub>-124.0</sub></td><td>1.77<sup>+1.1</sup><sub>-1.7</sub></td><td>36.7<sup>+22.9</sup><sub>-32.7</sub></td></tr><tr><td>YMR<sup>a</sup></td><td>29.8<sup>+18.9</sup><sub>-20.1</sub></td><td>13.5<sup>+6.4</sup><sub>-5.8</sub></td><td>80.8<sup>+6.9</sup><sub>-14.3</sub></td></tr><tr><td>OMP<sup>a</sup></td><td>-40.1<sup>+78.5</sup><sub>-134.1</sub></td><td>2.0<sup>+1.4</sup><sub>-1.8</sub></td><td>29.8<sup>+19.5</sup><sub>-33.2</sub></td></tr></tbody></table>

**Note.**

The N-body system was modeled using a total number of $N = 50,000$ ($N_{\text{OMP}} = 41,786$ and $N_{\text{YMR}} = 8214$) of single-mass particles. The mass of each particle is $m_* \approx 28 M_{\odot}$, and we adopted a softening length of 0.01 pc to smooth the close encounters between particles. Taking into account both of these approximations, the simulation time was rescaled to the evolutionary time of the system formed by the actual number of stars, assuming an overall average stellar mass of 0.94 $M_{\odot}$ (see Section 4.4), by using the ratio between the relaxation times of the real and simulated system as described in Mastrobuono-Battisti & Perets (2013). The initial conditions were built using the NEMO toolkit (Teuben 1995). We evolved the system for 2 Gyr—the age of the YMR population—using a version of the direct N-body code phiGRAPE adapted to run on GPUs (Harfst et al. 2007; Gaburov et al. 2009). The coevolution of the two populations and, in particular, the relaxation of the initially flattened YMR component, modify the final shape and kinematics of the whole cluster.

To evaluate the ellipticity of the different populations at the end of the simulation, we calculated their axial ratios using the inertia moments, as detailed in Katz (1991). The $b/a$ ratio between the intermediate and major axis is approximately equal

<footnote><sup>a</sup> For comparison, we add the values obtained assuming a distance to M54 of 28.4 kpc (Siegel et al. 2011).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Result of the fit for the YMR stars shown for the three components of the velocity vector, namely the LOS component, the proper motion component parallel to the rotation axis, and the proper motion component perpendicular to the rotation axis. The red line indicates the result of the MCMC, and the gray lines are 100 samples from the chain.</figcaption></figure>

to 1.0 at all radii, and the ellipticity is defined as $\epsilon = 1 - c/a$, where c is the minor axis of the system. The final ellipticity of the YMR population varies significantly with the distance from the center of the cluster; after 2 Gyr of evolution, the YMR stars are still in a significantly flattened configuration with ellipticity close to the initial one. In particular, the young component is almost spherical in the central regions of the cluster, while the ellipticity increases up to $\epsilon = 0.55$ at radii larger than 5 pc. The ellipticity is $\approx 0.4$ at 7 pc, which is the half-mass–radius of the cluster. The slight decrease in the ellipticity of the young population corresponds to an increase of flattening in the old stellar component. After 2 Gyr of coevolution with the younger component, the initially spherical OMP population becomes more flattened, reaching a maximum ellipticity of around 0.1. The ellipticity of this population decreases steadily with the distance from the cluster center (see left panel of Figure 8).

The change in morphology corresponds to a slight increase in the velocity anisotropy, whose amount is parameterized by the quantity $\beta_z = 1 - (\sigma_z/\sigma_R)^2$ in cylindrical coordinates with the z axis parallel to the angular momentum vector of the system (see right panel of Figure 8). Thus, $\beta_z$ is related to how the system is flattened (in this case the system is oblate, given the positive value of $\beta_z$), depending on the ratio between the dispersion in the z direction and in the plane perpendicular to the total angular momentum of the system. While, after 2 Gyr, the YMR disk still rotates with a peak velocity of around $11 \text{ km s}^{-1}$, the OMP population only shows a weak rotation pattern, with a peak velocity of $\sim 1 \text{ km s}^{-1}$ (see Figure 9). As already found by Mastrobuono-Battisti & Perets (2013, 2016), this result suggests that the coevolution of an initially spherical system with an embedded disk leads to an angular momentum redistribution between the flattened YMR and the OMP population. From this process follows the mixing of the

YMR and the OMP populations and the increased flattening of the OMP population (see left panel of Figure 10 for the angular momentum evolution of the two populations). The growing flattening is accompanied by only a slight increase in the rotation speed of the OMP population. Finally, the two populations, after 2 Gyr of evolution, are not yet fully mixed, with the YMR population still more centrally concentrated and less extended than the OMP population (see right panel of Figure 10).

## 7. Discussion

The kinematic results we present in this paper coupled to the dynamical simulations, add substantial evidence to the scenario we propose in Paper I from the findings based on the stellar characterization analysis on M54 and its populations.

### 7.1. Kinematic Comparison

As we mention in Paper I and show in detail in this work, all populations are spatially coincident, as evidenced by their measured spatial and LOS velocity distributions (see also Da Costa & Armandroff 1995; Ibata et al. 1997; Monaco et al. 2005; Bellazzini et al. 2008). We observe that the three populations display different velocity dispersion profiles. For the OMP population, the observed velocity dispersion is well described by a Plummer profile, with a median central velocity dispersion of $\sigma_0 = 15.30 \pm 0.54 \text{ km s}^{-1}$ and $\sigma \sim 10 \text{ km s}^{-1}$ at $r = 1'6$. The profiles are close to flat for the YMR and IMR populations. These are consistent with the velocity dispersion estimates from previous studies based on hundreds of stars (e.g., Bellazzini et al. 2008; Baumgardt & Hilker 2018). The flat velocity dispersion profile observed for the YMR population could be explained by (i) energy equipartition (see Section 4.4); (ii) the formation of the YMR population in a dynamically cold component, which is consistent with being rotationally supported and with low-velocity dispersion; (iii) crowding effects due to the highly crowded central region of M54 and the dominance of the OMP population over the YMR at these observed wavelengths. These effects limit our ability to extract an accurate number of stars that can produce an apparent drop in the central velocity dispersion. We are working on chemo-dynamical models on M54 that will help us to better constrain these kinematic differences between populations considering, e.g., their different spatial distribution and amount of rotation.

The IMR population displays a central velocity dispersion consistent with the OMP population; however, the profile is rather flat at a value of $\sim 15 \text{ km s}^{-1}$. This value at the largest radii is slightly higher than that found between $2' < r < 8'$ by Bellazzini et al. (2008), who showed the Sgr dSph stars followed a flat velocity dispersion profile at $\sim 10 \text{ km s}^{-1}$. This difference may be due to our particular definition of the IMR population, which excludes younger stars (at preferentially low dispersion at large radii) that may be in the Bellazzini et al. (2008) sample. At larger radii until $\sim 60'$, the Sgr dSph velocity dispersion profile (see Majewski et al. 2013) seems to be consistent with other dSph galaxies around the Milky Way (e.g., Walker et al. 2007). The flat dispersion profile of the IMR and in the outskirts of Sgr dSph, and its similarity to other dSphs could have several origins. For example in classical dSphs, this is interpreted as a result of large dark-matter-dominated potentials; alternatively, it could be due to tidally
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Ellipticity (left panel) and anisotropy parameter $\beta_z$ (right panel) as a function of radius for the whole system (gray), the OMP (blue), and YMR (red) populations at 0 Gyr (dashed lines) and after 2 Gyr (solid lines) of evolution as obtained from the *N*-body simulations. The angular momentum lost by the YMR population is acquired by the OMP population, that becomes slightly flattened and acquires a small amount of velocity anisotropy along the *z* axis. Plots are obtained considering the edge-on view of the cluster.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Rotation curves for the whole system (gray), the OMP (blue), and YMR (red) populations at 0 Gyr (dashed lines) and after 2 Gyr (solid lines) of evolution as obtained from the *N*-body simulations. While the YMR population still rotates significantly, the OMP population has acquired a small (~1 km s⁻¹) rotational speed. The cluster is seen edge-on.</figcaption></figure>

stripped or potential escaper stars (Bianchini et al. 2019). In such a complex system, the details of tidal disruption, including how much baryonic and dark matter remain in the various volumes and the joint effect on velocity dispersion profiles, are far from certain for this specific case. A large radial coverage of Sgr dSph along with detailed simulations is likely needed to understand these processes. Unfortunately, we cannot constrain these effects since our data is confined to the innermost 2', corresponding to the nuclear region, which might or might not be affected by tidal effects yet.

Our large MUSE data set allowed us to extract precise kinematics by using a large number of stars for each population. We find that the YMR population rotates ($A_{\text{rot}} = 4.91 \pm 0.90 \text{ km s}^{-1}$ and $V_{\text{max}} = 7.98 \pm 1.01 \text{ km s}^{-1}$) at a considerably higher speed than the OMP, which shows a weak sign of rotation ($A_{\text{rot}} = 0.75 \pm 0.55 \text{ km s}^{-1}$ and $V_{\text{max}} = 1.58 \pm 1.14 \text{ km s}^{-1}$). Note that the rotation obtained from 3D kinematics (see Table 3) is higher than the one obtained with the LOS data only (yet, consistent within the uncertainties). This difference might be the result of the subsample of stars from the MUSE sample for the estimate using the *Gaia* data where we obtain large error bars. Additionally, the 3D case does not find exactly an edge-on angle (~80°), which can increase the intrinsic $A_{\text{rot}}$.

We note that *Gaia* proper motions are affected by systematic errors, which have the same order of magnitude as the rotation signal we are trying to measure (0.07 mas yr⁻¹, 8–9 km s⁻¹ at M54’s distance). For M54, the uncertainties of *Gaia* proper motions are fairly large, due to the rather large distance and the high stellar density at the center. This reduces the accuracy of proper motions measurements, which are typically found to have errors of the order of ~20 km s⁻¹. All of these factors can effectively hide any rotation signal present in the plane of the sky. Future *Gaia* data releases are needed to obtain more accurate 3D rotation results. There are *Hubble Space Telescope* proper motions for M54 (Bellini et al. 2014; Watkins et al. 2015); however, they are relative, not absolute measurements, and cannot be used to study rotation.

Using *Gaia DR2* data of stars in both populations, we found that the OMP and YMR populations are comoving in 3D space. This suggests these two populations are dynamically bound, adding substantial evidence for spatial coincidence and against chance alignment.

## 7.2. YMR: Evidence of In Situ Formation

Young sub-components in NSCs have been detected in different types of galaxies (e.g., early- and late-type). These young structures are more centrally concentrated and flattened
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Left panel: evolution of the average angular momentum per unit mass perpendicular to the maximum rotation plane for the YMR (red) and OMP (blue) populations. The angular momentum lost by the YMR stars is redistributed among the OMP stars that, consequently, acquire a coherent rotational pattern and settle on a slightly flattened configuration. Right panel: spatial density profiles of the system, from the *N*-body simulations, considered as a whole (gray) and of the YMR (red) and OMP (blue) components taken separately as obtained from the *N*-body simulation.</figcaption></figure>

(Seth et al. 2006; Carson et al. 2015; Feldmeier-Krause et al. 2015; Nguyen et al. 2017). A similar central concentration was found in ω Cen (van de Ven et al. 2006), a presumed former NSC.

In Paper I, we found evidence to suggest that the YMR population formed ~3 Gyr ago likely in situ in a star-forming gas disk. In addition, we found that this population is the most centrally concentrated and the most flattened ($\epsilon = 0.31 \pm 0.10$).

In this paper, we add that the YMR population has a rotation amplitude of ~5 km s⁻¹, reaching a maximum rotation of ~8 km s⁻¹, consistent with the high degree of flattening that we found. In contrast, the OMP population shows a small rotation of 0.7 km s⁻¹, but it shows an amount of flattening similar to other clusters with higher rotation ($\epsilon \sim 0.15$). We will discuss this in the next section. From *N*-body simulations, Mastrobuono-Battisti & Perets (2013, 2016) and Mastrobuono-Battisti et al. (2019b) found that the young structures can survive despite being embedded in an older population. Moreover, they are found to rotate faster than the population in which they are embedded. Given the significantly different kinematics of these populations, it appears likely that the small difference observed in ellipticity is real. A representative case of a young stellar disk is observed in the Milky Way NSC (Genzel et al. 2003; Paumard et al. 2006; Bartko et al. 2009; Lu et al. 2009, 2013; Yelda et al. 2014; Feldmeier-Krause et al. 2015).

### 7.3. Kinematic Effects on the OMP Stars by the YMR

The formation and kinematics of the YMR population can affect the kinematics and shape of the population in which they were born. We found that the OMP population has a low rotation of $A_{\text{rot}} = 0.75 \pm 0.55$ km s⁻¹, with an ellipticity of $\epsilon = 0.16 \pm 0.06$. In the $(V/\sigma_0, \epsilon)$ diagram (see Figure 5), we observe that the OMP population (blue square) is more flattened than most of the GCs rotating at a similar speed and is similar to the most flattened clusters, e.g., ω Cen (black

diamond), M22 (black pentagon), which seem to rotate even faster.

Our *N*-body simulations consider a disk-like centrally concentrated YMR population, embedded in an initially spherical and five-times-more-massive OMP population. The two populations coevolve for 2 Gyr, which is the estimated age of the young population observed in M54. At the end of the simulation, the YMR population, which is born in a disk because of the angular momentum of the progenitor gas, relaxes and evolves toward a more spherical configuration while redistributing its angular momentum among the OMP stars. As a consequence of this process, the YMR disk slows down its rotation. At the same time, the OMP population acquires the angular momentum lost by the younger stars, decreases its velocity dispersion in the z direction and becomes slightly more flattened. This phenomenon could, at least partially, explain the ellipticity observed for the OMP population. Previous stellar populations, born in a disk, could have also contributed to the flattening of the OMP population leading to a more significant final effect, and shaping the present morphology of the OMP population. The relative angular momentum and ellipticity evolution due to the YMR interacting with the OMP population might be even stronger than our simulations suggest, given that in Section 4.4, the average mass of stars in the YMR population is approximately two times more than that in the OMP population, compared to the 1:1 ratio adopted in our simulations. Further simulations will help to address this effects.

Past mergers could have increased the ellipticity of the system but would have also significantly affected its kinematics, producing a higher rotational signal than that observed in M54 (Tsatsi et al. 2017; Mastrobuono-Battisti et al. 2019a). If the YMR population formed in a central disk that relaxed, leading to the observed OMP population flattening, we predict a radially increasing velocity anisotropy for the YMR population, with $\beta_z$ down to values smaller than 0.6 inside the half-mass-radius of the cluster (see Figure 8). The OMP population has a slight radial increase in the anisotropy profile.
Future Gaia observations will be able to verify these predictions at least in the outskirts of the cluster.

In Paper I, we found different spatial distributions for the OMP and YMR populations, with the YMR being the most centrally concentrated. Mastrobuono-Battisti & Perets (2013, 2016) also found that even after 12 Gyr, the two structures are not yet fully mixed, i.e., are still showing different distributions. In the simulations presented in this work, we also observe that after 2 Gyr of evolution (close to the age of the YMR population), the two populations are not yet fully mixed, and the YMR population is found to be still more centrally concentrated than the OMP population.

Note that the N-body simulations are not meant to exactly reproduce the current properties of M54 but rather to test the possibility that the interactions between the OMP and YMR populations might have mutually modified the dynamical configuration of these two populations.

### 7.4. OMP: Remnant of a Cluster Merger?

In Paper I, we found that the OMP population displays a large spread in both metallicity and age. A large spread in metallicity alone can be explained by self-enrichment during formation (Bailin 2018), but it does not explain the large spread in age. Thus, we suggested that this could be the result of a merger event between two or more clusters that fell into the central region of the host, as suggested through simulations (e.g., Amaro-Seoane et al. 2013; Bekki & Tsujimoto 2016; Gavagnin et al. 2016; Khoperskov et al. 2018; Mastrobuono-Battisti et al. 2019a). In addition, a contribution to this high spread in metallicity might be due to the old, metal-poor population of Sgr dSph that mixed with the merged GCs to the resulting OMP population of this NSC.

According to the N-body simulation presented in this work, the formation and evolution of the YMR population alone can explain the shape and kinematics of the OMP. However, it does not explain the large spread in both age and metallicity, that we found in this population.

The final rotation of two merged clusters strongly depends on the conditions, as the orbital configuration and relaxation states of the merging clusters do not always result in a highly rotating structure (Mastrobuono-Battisti et al. 2019a). Based on the kinematics, we cannot be certain that two merging clusters were actually involved in the formation of the OMP population but cannot discard this possibility either. Further high-precision observations and more detailed simulations appropriate to mergers in the Sgr potential may help to shed light on whether this process may have contributed to the build-up of the OMP population in the Sgr dSph.

### 7.5. Comparison with ω Cen

The similarities between NSCs in dwarf galaxies and high-mass, metal complex GCs suggest that such GCs might be former nuclei of dwarf galaxies accreted by the Milky Way (Zinnecker et al. 1988; Böker 2008; Da Costa 2016). One of the most likely former nuclei is ω Cen, which is the most massive GC in the Milky Way ((3.55 ± 0.03) × 10⁶ M☉, Baumgardt & Hilker 2018). It has a large spread in both iron (Johnson & Pilachowski 2010) and age (Hilker et al. 2004; Villanova et al. 2014) among its stars and hosts multiple stellar populations (Bellini et al. 2017; Milone et al. 2017), which present different kinematics (Bellini et al. 2018). These

properties, and their similarity to what we have found in M54, make a strong case for ω Cen to be a stripped NSC that once resided in a dwarf galaxy now accreted by the Milky Way (Lee et al. 1999; Majewski et al. 2000, 2012; Bekki & Freeman 2003; Carretta et al. 2010b). Using N-body simulations, Ibata et al. (2019a) found that the “Fimbulthul” structure detected with the Gaia DR2 observations (Ibata et al. 2019b) could be a tidal stream of ω Cen. This suggests that there may be stars in the halo that were once part of ω Cen; finding these escaped stars would help us to understand the progenitor. Massari et al. (2019) used kinematic information of the Milky Way GCs from Gaia to suggest that ω Cen could be the former NSC of the progenitor galaxy of the merger event Gaia-Enceladus (Helmi et al. 2018). Observations of ω Cen and M54 suggest that ω Cen could be in a more advanced “stripped nuclei stage” (see Pfeffer & Baumgardt 2013) than M54 (e.g., Bellazzini et al. 2008; Carretta et al. 2010b), since M54 is still observed at the photometric center of its host galaxy (Ibata et al. 1994), which is in ongoing disruption by the tidal field of the Milky Way (Ibata et al. 1997).

Similar to other GCs, ω Cen presents ellipticity variations as a function of radius (Geyer et al. 1983) and is found to be one of the most flattened Milky Way GCs (Meylan 1987). As we mentioned in Section 4.3, there is a correlation between the amount of rotation of a cluster and its ellipticity (Fabricius et al. 2014; Kamann et al. 2018), which has been detected for ω Cen (e.g., Meylan & Mayor 1986; Bianchini et al. 2013, 2018; Kamann et al. 2018). ω Cen behaves close to an isotropic oblate rotator in the inner parts ($< 10^{10}$), becoming more anisotropic at larger radii, probably due to the tidal effects by the Milky Way (van de Ven et al. 2006), which do not affect the inner parts of the cluster. For M54, we find a difference of $Δε ≈ 0.05$ between the ellipticity estimated over the entire FOV (~2.5 × $R_{HL}$) and at the half-light radius. Being cautious with the ellipticity uncertainties, this suggests that the ellipticity for M54 varies as a function of radius, becoming more flattened at larger radii. From our ($V/σ_0, ε$) diagram (see Figure 5), we observe that M54 is close to the isotropic oblate rotator relation when looking at its half-light radius, showing low rotation with a small degree of flattening. Watkins et al. (2015) used proper motions to show that M54 is isotropic out to the half-light radius. However, this does not imply that the cluster cannot display anisotropy at larger radii. For the observations of M54, further analysis and information are needed to confirm and constrain its presumed radial anisotropy.

Additionally, a disk-like component was detected in ω Cen by van de Ven et al. (2006). We found a similar structure in M54 corresponding to the YMR population that shows a strong rotational signal and that most likely formed in the center of the OMP population. We also speculate that the disk-like component in ω Cen was born within the old population, which can be tested by measuring its metallicity, age, and rotational structure.

Although small differences exist between ω Cen and M54, the similarity of their stellar populations, their morphology, and kinematics provides additional evidence that ω Cen is, in fact, the nucleus of a stripped dwarf galaxy. This comparison considers M54 as the NSC of Sgr dSph, which is now composed of three main populations. This nuclear stellar structure is what we would observe if its host galaxy, the Sgr dSph, was totally stripped.
## 8. Conclusions

In Paper I, we identified three subpopulations in M54, the NSC of the Sgr dSph based only on their ages and metallicities. They are: (i) YMR: 2.2 Gyr and [Fe/H] = -0.04, (ii) IMR: 4.3 Gyr and [Fe/H] = -0.29, and (iii) OMP: 12.2 Gyr and [Fe/H] = -1.41. Here, we have presented a kinematic analysis of each subpopulation to both characterize the populations and gain new insight into their formation and evolution.

In this paper, we continued the characterization of the stellar populations in M54 by adding kinematic information. M54 is one of the closest extragalactic NSC, at a distance where stars can be resolved and characterized, enabling the distinction of the stellar populations. In studies of extragalactic NSCs, the different stellar populations cannot be easily resolved, restricting the kinematic analysis to the nucleus as one single structure. In M54, we have the opportunity to directly compare the kinematics for the nucleus considering all of the stars in the sample (i.e., like one single structure), as is usually done for extragalactic NSCs, and for the different populations. We present the kinematics for two cases: (i) all M54 star members with no distinction, and (ii) for the populations detected in M54 in Paper I.

We found that all three populations show different velocity dispersion profiles; the YMR and IMR populations follow a flat distribution, whereas the OMP population is well-defined by a Plummer profile. We find that all populations show a different amount of rotation. The YMR population shows a considerable amount of rotation (~5 km s⁻¹), followed by the IMR population, which also rotates but more slowly (~2 km s⁻¹). We detect a weak signal of rotation for the OMP population (<1 km s⁻¹). Hence, the findings from the kinematic analysis support our findings from Paper I that the populations do not have a common origin.

From these findings together with those in Paper I, we suggest the following conclusions:

1. From our large sample of stars, we find that all populations have the same systemic velocity. Combining our data with proper motions from Gaia DR2, we find that the stars in the OMP and YMR population are comoving in 3D space. Taken together, this secures the finding that the populations are spatially coincident and disfavors chance alignment.

2. The YMR population displays a high amount of rotation. Combined with the fact that it is more flattened and more centrally concentrated, this strongly favors the scenario wherein these stars formed in situ from enriched gas accreted at the center of M54.

3. The OMP population appears more flattened than typical GCs rotating at the same speed. The N-body simulation emulating the YMR-OMP system in M54 suggests that this could result from the angular momentum transferred from the YMR to the OMP population, thus decreasing its velocity dispersion in the z direction and becoming more flattened. However, this alone does not explain the high iron and age spread this population presents, which we suggested to be signs of a merger event between two or more GCs in the center of the host. The merger scenario cannot be ruled out since mergers of GCs both do or do not end in a rotating structure depending on the initial conditions. The OMP population might also be the result of a GC that sank to the center where the old, metal-poor population from Sgr dSph was already in place, thus together contributing to the high spread in age and metallicity.

4. We tested the possibility that effects connected to energy equipartition could imprint different kinematics due to energy exchange between stars with different masses. We find that the difference in velocity dispersions between the YMR and OMP populations are consistent with equipartition effects driven by relaxation processes.

5. The current information we obtain from the IMR population does not allow us to constrain if its stars are dynamically bound to the nucleus. However, the kinematics of this population are consistent with these stars being part of the central field of the Sgr dSph.

The kinematic analysis presented in this paper certainly adds an essential piece of information to the understanding of M54 and its populations.

NSCs display multiple populations—as observed in the nucleus of the Milky Way—and grow with time depending on their environmental conditions. The nucleus of the Sgr dSph possibly started with an old and metal-poor population where one or more GCs merged, thus forming together an old, metal-poor stellar structure with a considerable spread in age and metallicity. In spite of our limited ability to differentiate them, this is what we observe today as the OMP component of the nucleus of the Sgr dSph. Then, ~3 Gyr ago a young population (YMR) was born within the nucleus. All of these populations in this NSC that Messier observed for the first time in 1778 and named M54 form the NSC in the Sgr dSph. This structure with all its populations is likely what we would observe after the total stripping of the host, as in the case of ω Cen.

## 9. Future Prospects

The confirmation of the putative IMBH at the center of M54 would provide a strong piece of evidence to support its nuclear origin. It would also imply their existence in other NSCs in low-mass galaxies where they have not been detected, possibly due to instrumentation limitations. Many high-mass metal complex GCs are IMBH host candidates, but there is no confirmation yet, e.g., ω Cen (Noyola et al. 2008; Jalali et al. 2012; Baumgardt 2017), G1 (Gebhardt et al. 2005), NGC 6388 (Lützgendorf et al. 2011). M54 is presumed to host an IMBH of $10^4 M_{\odot}$ (Ibata et al. 2009; Baumgardt 2017), but a confirmation is still needed.

We are performing dynamical modeling of M54 and its individual stellar populations to constrain the existence and mass of the IMBH in its center. This will also constrain variations in anisotropy as a function of radius. To achieve this, we work with an additional MUSE narrow-field mode data set that allows us to extract hundreds of stars in the innermost 3" of M54, in the sphere of influence of the presumed IMBH. In an upcoming paper, we will present discrete Jeans dynamical models for M54 where we expect to better constrain its dynamical properties and the existence of the IMBH.

S.K. gratefully acknowledges funding from a European Research Council consolidator grant (ERC-CoG-646928- Multi-Pop). R.L. acknowledges funding from a Natural Sciences and Engineering Research Council of Canada PDF award, and A.M.B., R.L., and N.N. acknowledge support by
Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject A7, A8 and B8)—Project-ID 138713538—of the German Research Foundation (DFG). M.A.C., R.L., and P.B. acknowledge support from DAAD PPP project number 57316058 “Finding and exploiting accreted star clusters in the Milky Way.” Work by A.C.S. on this project was supported by NSF grant AST-1350389. G.v.d.V. and L.L.W. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 724857 (Consolidator Grant ArcheoDyn). This research has made use of NASA’s Astrophysics Data System. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

## ORCID iDs

M. Alfaro-Cuello https://orcid.org/0000-0002-1212-2844

N. Kacharov https://orcid.org/0000-0002-6072-6669

N. Neumayer https://orcid.org/0000-0002-6922-2598

P. Bianchini https://orcid.org/0000-0002-0358-4502

A. Mastrobuono-Battisti https://orcid.org/0000-0002-2386-9142

A. C. Seth https://orcid.org/0000-0003-0248-5470

T. Böker https://orcid.org/0000-0002-5666-7782

S. Kamann https://orcid.org/0000-0001-6604-0505

L. L. Watkins https://orcid.org/0000-0002-1343-134X

G. van de Ven https://orcid.org/0000-0003-4546-7731

## References

Alfaro-Cuello, M., Kacharov, N., Neumayer, N., et al. 2019, ApJ, 886, 57

Amaro-Seoane, P., Konstantinidis, S., Brem, P., & Catelan, M. 2013, MNRAS, 435, 809

Bacon, R., Vernet, J., Borisova, E., et al. 2014, Msngr, 157, 13

Bailin, J. 2018, ApJ, 863, 99

Bartko, H., Martins, F., Fritz, T. K., et al. 2009, ApJ, 697, 1741

Baumgardt, H. 2017, MNRAS, 464, 2174

Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520

Bekki, K., & Freeman, K. C. 2003, MNRAS, 346, L11

Bekki, K., & Tsujimoto, T. 2016, ApJ, 831, 70

Bellazzini, M., Bragaglia, A., Carretta, E., et al. 2012, A&A, 538, A18

Bellazzini, M., Ibata, R. A., Chapman, S. C., et al. 2008, AJ, 136, 1147

Bellini, A., Anderson, J., van der Marel, R. P., et al. 2014, ApJ, 797, 115

Bellini, A., Libralato, M., Bedin, L. R., et al. 2018, ApJ, 853, 86

Bellini, A., Milone, A. P., Anderson, J., et al. 2017, ApJ, 844, 164

Bianchini, P., Ibata, R., & Famaey, B. 2019, ApJL, 887, L12

Bianchini, P., van de Ven, G., Norris, M. A., Schinnerer, E., & Varri, A. L. 2016, MNRAS, 458, 3644

Bianchini, P., van der Marel, R. P., del Pino, A., et al. 2018, MNRAS, 481, 2125

Bianchini, P., Varri, A. L., Bertin, G., & Zocchi, A. 2013, ApJ, 772, 67

Binney, J. 2005, MNRAS, 363, 937

Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, NJ: Princeton Univ. Press)

Böker, T. 2008, ApJL, 672, L111

Böker, T., Laine, S., van der Marel, R. P., et al. 2002, AJ, 123, 1389

Böker, T., Sarzi, M., McLaughlin, D. E., et al. 2004, AJ, 127, 105

Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345

Cappellari, M., Emsellem, E., Bacon, R., et al. 2007, MNRAS, 379, 418

Carollo, C. M., Stiavelli, M., & Mack, J. 1998, AJ, 116, 68

Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2010a, A&A, 520, A95

Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2010b, ApJL, 714, L7

Carson, D. J., Barth, A. J., Seth, A. C., et al. 2015, AJ, 149, 170

Cordero, M. J., Hénault-Brunet, V., Pilachowski, C. A., et al. 2017, MNRAS, 465, 3515

Côté, P., Piatek, S., Ferrarese, L., et al. 2006, ApJS, 165, 57

Da Costa, G. S. 2016, in IAU Symp. 317, The General Assembly of Galaxy Halos: Structure, Origin and Evolution, ed. A. Bragaglia et al. (Cambridge: Cambridge Univ. Press), 110

Da Costa, G. S., & Armandroff, T. E. 1995, AJ, 109, 2533

den Brok, M., Peletier, R. F., Seth, A., et al. 2014, MNRAS, 445, 2385

Einsel, C., & Spurzem, R. 1999, MNRAS, 302, 81

Emsellem, E., Cappellari, M., Krajnović, D., et al. 2011, MNRAS, 414, 888

Fabricius, M. H., Noyola, E., Rukdee, S., et al. 2014, ApJL, 787, L26

Feldmeier, A., Neumayer, N., Seth, A., et al. 2014, A&A, 570, A2

Feldmeier-Krause, A., Neumayer, N., Schödel, R., et al. 2015, A&A, 584, A2

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaburov, E., Harfst, S., & Portegies Zwart, S. 2009, NewA, 14, 630

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Gavagnin, E., Mapelli, M., & Lake, G. 2016, MNRAS, 461, 1276

Gebhardt, K., Rich, R. M., & Ho, L. C. 2005, ApJ, 634, 1093

Genzel, R., Schödel, R., Ott, T., et al. 2003, ApJ, 594, 812

Georgiev, I. Y., & Böker, T. 2014, MNRAS, 441, 3570

Geyer, E. H., Hopp, U., & Nelles, B. 1983, A&A, 125, 359

Gibbons, S. L. J., Belokurov, V., & Evans, N. W. 2017, MNRAS, 464, 794

Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65

Harfst, S., Gualandris, A., Merritt, D., et al. 2007, NewA, 12, 357

Harris, W. E. 1996, AJ, 112, 1487

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85

Hénault-Brunet, V., Gieles, M., Agertz, O., & Read, J. I. 2015, MNRAS, 450, 1164

Hilker, M., Kayser, A., Richtler, T., & Willemsen, P. 2004, A&A, 422, L9

Ibata, R., Bellazzini, M., Chapman, S. C., et al. 2009, ApJL, 699, L169

Ibata, R., Bellazzini, M., Malhan, K., Martin, N., & Bianchini, P. 2019a, NatAs, 3, 667

Ibata, R., Gilmore, G., & Irwin, M. J. 1994, Natur, 370, 194

Ibata, R. A., Malhan, K., & Martin, N. F. 2019b, ApJ, 872, 152

Ibata, R. A., Wyse, R. F. G., Gilmore, G., Irwin, M. J., & Suntzeff, N. B. 1997, AJ, 113, 634

Jalali, B., Baumgardt, H., Kissler-Patig, M., et al. 2012, A&A, 538, A19

Jindal, A., Webb, J. J., & Bovy, J. 2019, MNRAS, 487, 3693

Johnson, C. I., & Pilachowski, C. A. 2010, ApJ, 722, 1373

Kacharov, N., Bianchini, P., Koch, A., et al. 2014, A&A, 567, A69

Kacharov, N., Neumayer, N., Seth, A. C., et al. 2018, MNRAS, 480, 1973

Kamann, S., Husser, T.-O., Dreizler, S., et al. 2018, MNRAS, 473, 5591

Katz, N. 1991, ApJ, 368, 325

Khoperskov, S., Mastrobuono-Battisti, A., Di Matteo, P., & Haywood, M. 2018, A&A, 620, A154

Kimmig, B., Seth, A., Ivans, I. I., et al. 2015, AJ, 149, 53

King, I. R. 1966, AJ, 71, 64

Koch, A., Hanke, M., & Kacharov, N. 2018, A&A, 616, A74

Koleva, M., Prugniel, P., Bouchard, A., & Wu, Y. 2009, A&A, 501, 1269

Kruijssen, J. M. D., Pfeffer, J. L., Reina-Campos, M., Crain, R. A., & Bastian, N. 2019, MNRAS, 486, 3180

Lardo, C., Pancino, E., Bellazzini, M., et al. 2015, A&A, 573, A115

Lee, Y.-W., Joo, J.-M., Sohn, Y.-J., et al. 1999, Natur, 402, 55

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Lu, J. R., Do, T., Ghez, A. M., et al. 2013, ApJ, 764, 155

Lu, J. R., Ghez, A. M., Hornstein, S. D., et al. 2009, ApJ, 690, 1463

Lützgendorf, N., Kissler-Patig, M., Noyola, E., et al. 2011, A&A, 533, A36

Mackey, A. D., Da Costa, G. S., Ferguson, A. M. N., & Yong, D. 2013, ApJ, 762, 65

Majewski, S. R., Hasselquist, S., Łokas, E. L., et al. 2013, ApJL, 777, L13

Majewski, S. R., Nidever, D. L., Smith, V. V., et al. 2012, ApJL, 747, L37

Majewski, S. R., Patterson, R. J., Dinescu, D. I., et al. 2000, in Liege Int. Astrophysical Colloquia 35, The Galactic Halo: From Globular Cluster to Field Stars, ed. A. Noels et al. (Liege: Institut d'Astrophysique et de Geophysique), 619

Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4

Mastrobuono-Battisti, A., Khoperskov, S., Di Matteo, P., & Haywood, M. 2019a, A&A, 622, A86

Mastrobuono-Battisti, A., & Perets, H. B. 2013, ApJ, 779, 85

Mastrobuono-Battisti, A., & Perets, H. B. 2016, ApJ, 823, 61

Mastrobuono-Battisti, A., Perets, H. B., Gualand ris, A., Neumayer, N., & Sippel, A. C. 2019b, MNRAS, 490, 5820

Mylan, G. 1987, A&A, 184, 144
Meylan, G., & Mayor, M. 1986, A&A, **166**, 122

Milone, A. P., Marino, A. F., Bedin, L. R., et al. 2017, MNRAS, **469**, 800

Misgeld, I., & Hilker, M. 2011, MNRAS, **414**, 3699

Monaco, L., Bellazzini, M., & Bonifacio, P. 2005, A&A, **441**, 141

Mucciarelli, A., Bellazzini, M., Ibata, R., et al. 2017, A&A, **605**, A46

Muñoz, R. P., Eigenthaler, P., Puzia, T. H., et al. 2015, ApJL, **813**, L15

Neumayer, N., Seth, A., & Boeker, T. 2020, arXiv:2001.03626

Nguyen, D. D., Seth, A. C., den Brok, M., et al. 2017, ApJ, **836**, 237

Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2018, ApJ, **858**, 118

Niederste-Ostholt, M., Belokurov, V., & Evans, N. W. 2012, MNRAS, **422**, 207

Norris, M. A., Kannappan, S. J., Forbes, D. A., et al. 2014, MNRAS, **443**, 1151

Noyola, E., Gebhardt, K., & Bergmann, M. 2008, ApJ, **676**, 1008

Paumard, T., Genzel, R., Martins, F., et al. 2006, ApJ, **643**, 1011

Pfeffer, J., & Baumgardt, H. 2013, MNRAS, **433**, 1997

Pfeffer, J., Griffen, B. F., Baumgardt, H., & Hilker, M. 2014, MNRAS, **444**, 3670

Phillips, A. C., Illingworth, G. D., MacKenty, J. W., & Franx, M. 1996, AJ, **111**, 1566

Plummer, H. C. 1911, MNRAS, **71**, 460

Sánchez-Janssen, R., Côté, P., Ferrarese, L., et al. 2019, ApJ, **878**, 18

Seth, A. C., Dalcanton, J. J., Hodge, P. W., & Debattista, V. P. 2006, AJ, **132**, 2539

Seth, A. C., van den Bosch, R., Mieske, S., et al. 2014, Natur, **513**, 398

Siegel, M. H., Dotter, A., Majewski, S. R., et al. 2007, ApJ, **667**, L57

Siegel, M. H., Majewski, S. R., Law, D. R., et al. 2011, ApJ, **743**, 20

Sollima, A., Baumgardt, H., & Hilker, M. 2019, MNRAS, **485**, 1460

Sollima, A., Bellazzini, M., Smart, R. L., et al. 2009, MNRAS, **396**, 2183

Spitzer, L. J. 1969, ApJL, **158**, L139

Tepper-García, T., & Bland-Hawthorn, J. 2018, MNRAS, **478**, 5263

Teuben, P. 1995, in ASP Conf. Ser. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes (San Francisco, CA: ASP), 398

Tiongco, M. A., Vesperini, E., & Varri, A. L. 2018, MNRAS, **475**, L86

Trenti, M., & van der Marel, R. 2013, MNRAS, **435**, 3272

Tsatsi, A., Mastrobuono-Battisti, A., van de Ven, G., et al. 2017, MNRAS, **464**, 3720

Turner, M. L., Côté, P., Ferrarese, L., et al. 2012, ApJS, **203**, 5

van de Ven, G., van den Bosch, R. C. E., Verolme, E. K., & de Zeeuw, P. T. 2006, A&A, **445**, 513

Vasiliev, E. 2019, MNRAS, **484**, 2832

Villanova, S., Geisler, D., Gratton, R. G., & Cassisi, S. 2014, ApJ, **791**, 107

Voggel, K. T., Seth, A. C., Baumgardt, H., et al. 2019, ApJ, **871**, 159

Walcher, C. J., van der Marel, R. P., McLaughlin, D., et al. 2005, ApJ, **618**, 237

Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2007, ApJL, **667**, L53

Watkins, L. L., van der Marel, R. P., Bellini, A., & Anderson, J. 2015, ApJ, **803**, 29

Wu, Y., Singh, H. P., Prugniel, P., Gupta, R., & Koleva, M. 2011, A&A, **525**, A71

Yelda, S., Ghez, A. M., Lu, J. R., et al. 2014, ApJ, **783**, 131

Zinnecker, H., Keable, C. J., Dunlop, J. S., Cannon, R. D., & Griffiths, W. K. 1988, in IAU Symp. 126, The Harlow-Shapley Symp. Globular Cluster Systems in Galaxies, ed. J. E. Grindlay & A. G. D. Philip (Dordrecht: Kluwer), 603