# A Deep View into the Nucleus of the Sagittarius Dwarf Spheroidal Galaxy with MUSE. III. Discrete Multicomponent Population-dynamical Models Based on the Jeans Equations

Nikolay Kacharov, Mayte Alfaro-Cuello, Nadine Neumayer, Nora Lützgendorf, Laura L. Watkins, Alessandra Mastrobuono-Battisti, Sebastian Kamann, Glenn van de Ven, Anil C. Seth, Karina T. Voggel, et al.

► To cite this version:

Nikolay Kacharov, Mayte Alfaro-Cuello, Nadine Neumayer, Nora Lützgendorf, Laura L. Watkins, et al.. A Deep View into the Nucleus of the Sagittarius Dwarf Spheroidal Galaxy with MUSE. III. Discrete Multicomponent Population-dynamical Models Based on the Jeans Equations. The Astrophysical Journal, 2022, 939, 10.3847/1538-4357/ac9280. insu-03865900

HAL Id: insu-03865900

https://insu.hal.science/insu-03865900

Submitted on 23 Nov 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A Deep View into the Nucleus of the Sagittarius Dwarf Spheroidal Galaxy with MUSE.
## III. Discrete Multicomponent Population-dynamical Models Based on the Jeans Equations

Nikolay Kacharov¹, Mayte Alfaro-Cuello² ID, Nadine Neumayer³ ID, Nora Lützgendorf⁴ ID, Laura L. Watkins⁵ ID,
Alessandra Mastrobuono-Battisti⁶ ID, Sebastian Kamann⁷, Glenn van de Ven⁸ ID, Anil C. Seth⁹ ID, Karina T. Voggel¹⁰ ID,
Iskren Y. Georgiev³ ID, Ryan Leaman³,⁸, Paolo Bianchini¹⁰ ID, Torsten Böker⁴ ID, and Steffen Mieske¹¹

¹Leibniz Institute for Astrophysics, An der Sternwarte 16, D-14482 Potsdam, Germany; kacharov@aip.de

²Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

³Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

⁴European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore MD 21218, USA

⁵AURA for the European Space Agency, ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore MD 21218, USA

⁶GEPI, Observatoire de Paris, PSL Research University, CNRS, Place Jules Janssen, F-92190 Meudon, France

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

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

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

¹⁰Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France

¹¹European Southern Observatory, 3107 Alonso de Córdova, Vitacura, Santiago, Chile

Received 2022 July 11; revised 2022 September 12; accepted 2022 September 12; published 2022 November 14

## Abstract

We present comprehensive multicomponent dynamical models of M54 (NGC 6715), the nuclear star cluster of the Sagittarius (Sgr) dwarf galaxy, which is undergoing a tidal disruption in the Milky Way halo. Previous papers in this series used a large MUSE mosaic data set to identify multiple stellar populations in the system and study their kinematic differences. Here, we use Jeans-based dynamical models that fit the population properties (mean age and metallicity), spatial distributions, and kinematics simultaneously. They provide a solid physical explanation for our previous findings. Population-dynamical models deliver a comprehensive view of the whole system, and allow us to disentangle the different stellar populations. We explore their dynamical interplay and confirm our previous findings about the build-up of Sgr's nuclear cluster via contributions from globular cluster stars, Sgr inner field stars, and in situ star formation. We explore various parameterizations of the gravitational potential and show the importance of a radially varying mass-to-light ratio for the proper treatment of the mass profile. We find a total dynamical mass within M54's tidal radius (~75 pc) of $1.60 \pm 0.07 \times 10^6 M_{\odot}$ in excellent agreement with N-body simulations. Metal-poor globular cluster stars contribute about 65% of the total mass or $1.04 \pm 0.05 \times 10^6 M_{\odot}$. Metal-rich stars can be further divided into young and intermediate-age populations, which contribute $0.32 \pm 0.02 \times 10^6 M_{\odot}$ (20%) and $0.24 \pm 0.02 \times 10^6 M_{\odot}$ (15%), respectively. Our population-dynamical models successfully distinguish the different stellar populations in Sgr's nucleus because of their different spatial distributions, ages, metallicities, and kinematic features.

Unified Astronomy Thesaurus concepts: Stellar kinematics (1608); Globular star clusters (656); Galaxy dynamics (591); Ultracompact dwarf galaxies (1734); Stellar populations (1622); Galaxy nuclei (609)

Supporting material: machine-readable table

## 1. Introduction

Nuclear star clusters (NSC) are common among low- and intermediate-mass galaxies of all morphological types, including dwarf galaxies with stellar masses lower than $10^9 M_{\odot}$ (see the recent review by Neumayer et al. 2020). Their extreme stellar densities ($>10^6 M_{\odot} \text{pc}^{-3}$ Hopkins et al. 2010), possible coincidence with supermassive black holes (Filippenko & Ho 2003; Seth et al. 2008; Graham & Spitler 2009; Neumayer & Walcher 2012; Nguyen et al. 2019), and extended star formation histories (SFHs) with the coexistence of multiple stellar populations (Rossa et al. 2006; Walcher et al. 2006; Norris et al. 2015; Kacharov et al. 2018; Fahrion et al. 2021; Hannah et al. 2021)

make such objects very interesting from a population-dynamical point of view.

The Sagittarius (Sgr) dwarf spheroidal galaxy (dSph), discovered by Ibata et al. (1994), is the closest known example of a nucleated dwarf galaxy (Bellazzini et al. 2008). It is currently being tidally stripped by the Milky Way, but its nucleus—the complex, massive, and dense star cluster M54 is still largely intact (Bassino et al. 1994; Bekki et al. 2003; Pfeffer & Baumgardt 2013; Pfeffer et al. 2021).

This work continues the series of papers on the Sgr dSph NSC—M54, based on Very Large Telescope (VLT) Multi-Unit Spectroscopic Explorer (MUSE) observations by Alfaro-Cuello et al. (2019) and Alfaro-Cuello et al. (2020), hereafter Paper I and Paper II. In Paper I, we introduced our observational sample and classified M54's stars in three different stellar populations, based on measured metallicities and ages. We named these populations old metal-poor (OMP), intermediate-age metal-rich (IMR), and young metal-rich

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
(YMR). We found that the YMR population is the most centrally concentrated and the IMR is the most spatially extended. In Paper II, we explored their kinematic properties and showed that they also differ considerably.

Sgr's nuclear formation picture that emerged from these studies is in agreement with the most widely accepted scenarios for NSC formation, namely, that both proposed mechanisms—in situ star formation from accreted gas (Mihos & Hernquist 1994; Milosavljević 2004; Schinnerer et al. 2008; Bekki 2015) and mergers of globular clusters (GCs; Tremaine et al. 1975; Oh & Lin 2000; Lotz et al. 2001; Capuzzo-Dolcetta & Miocchi 2008a, 2008b; Antonini et al. 2012; Gnedin et al. 2014) operate simultaneously (Neumayer et al. 2011; den Brok et al. 2014b; Antonini et al. 2015; Cole et al. 2017; Fahrion et al. 2022).

This paper is the third in a series of papers wherein we focus on constraining the overall shape of the gravitational potential and explore the interplay between the various distinct M54 stellar populations in the context of NSC formation theories.

The OMP population old age and kinematic properties are typical for GCs. They show negligible rotation and a radially decreasing velocity dispersion profile. Its spatial distribution can be well described with a King profile (King 1962; Trager et al. 1995; McLaughlin & van der Marel 2005; Monaco et al. 2005). Its significant metallicity spread (0.24 dex, Paper I; 0.19 dex Carretta et al. 2010), compared to the typical GC metallicity spread of <0.1 dex (Carretta et al. 2009), suggests that the OMP stars might be the result of multiple merged GCs that have sunk to Sgr's center via dynamical friction. In addition, stars from the metal-poor (MP) end of the Sgr field distribution likely also coexist in the nuclear cluster, inflating M54's OMP metallicity spread.

YMR stars, on the other hand, exhibit a significantly higher rate of systemic rotation of ~5 km s⁻¹ and have a lower velocity dispersion than OMP stars. This is expected if they formed in situ within a high angular momentum disk structure from dynamically cold gas that was accreted in Sgr's central region. Their mean age of ~2 Gyr coincides with Sgr's first peri-galactic passage, which could have triggered the nuclear star formation episode (Tepper-García & Bland-Hawthorn 2018; Di Cintio et al. 2021).

The IMR population has a significant age and metallicity spread and its SFH and the recovered age-metallicity relation (Paper I) follows that of Sgr's field population closely (Layden & Sarajedini 2000). This population is also characterized by a radially flatter velocity dispersion profile with higher dispersion in the outer nuclear region, compared to OMP stars. It is much more extended than the YMR and OMP populations; however, it still forms a density cusp, making it distinct from the field stars' spatial distribution (Monaco et al. 2005).

Previous dynamical studies of M54 have looked into the possibility of detecting an intermediate-mass black hole in the center of the system (Ibata et al. 2009; Wrobel et al. 2011), but have not reached conclusive results. Although our MUSE data is suitable to explore this question with improved sensitivity, we leave it for a future publication. In addition, as an integral part of the dark-matter dominated dSph galaxy, Sgr's nucleus is also likely embedded in a dark matter halo, whose contribution to the gravitational potential starts to become important in the outer regions (>30 pc, e.g., Carlberg & Grillmair 2022).

In this work, we base our dynamical models on the Jeans (1915) hydrodynamical equations of stellar motions, which are often used today due to their simplicity and computational

<table><caption>**Table 1**
M54 Adopted Parameters</caption><thead><tr><th>Parameter</th><th>Value<sup>a</sup></th></tr></thead><tbody><tr><td>R.A.</td><td>283°76387</td></tr><tr><td>decl.</td><td>-30°479861</td></tr><tr><td>PM<sub>R.A.</sub></td><td>-2.683 mas yr<sup>-1</sup></td></tr><tr><td>PM<sub>decl.</sub></td><td>-1.385 mas yr<sup>-1</sup></td></tr><tr><td>Distance</td><td>26.5 kpc</td></tr><tr><td>Core radius (r<sub>c</sub>)</td><td>5"4</td></tr><tr><td>Half-light radius (r<sub>h</sub>)</td><td>49"2</td></tr><tr><td>Tidal radius<sup>b</sup> (r<sub>t</sub>)</td><td>590"</td></tr><tr><td>Abs. V-band mag.</td><td>-9.98 mag</td></tr><tr><td>Luminosity (L<sub>v</sub>)</td><td>0.85 × 10<sup>6</sup> L<sub>☉</sub></td></tr><tr><td>MUSE r<sub>max</sub><sup>c</sup></td><td>149"</td></tr><tr><td>Angular scale</td><td>0.13 pc arcsec<sup>-1</sup></td></tr></tbody></table>

**Notes.**

<sup>a</sup> All values are taken from the Harris catalog of GCs (Harris 1996, 2010 edition), besides the proper motions, which come from Vasiliev & Baumgardt (2021).

<sup>b</sup> The concept of a tidal radius is not really meaningful in nuclear clusters. In any case, we give our total mass estimates at this radius.

<sup>c</sup> Diagonal extent of the MUSE mosaic field of view.

efficiency (Cappellari 2008). In the past, direct applications of this approach often required binning the data to infer the higher velocity moments, which significantly limited our ability to model complex stellar systems. Watkins et al. (2013) presented a method to use the Jeans dynamical models with discrete likelihood functions, which takes full advantage of the kinematic measurements for individual stars with their respective uncertainties and allows for the simultaneous modeling of multiple stellar populations (Zhu et al. 2016; Kamann et al. 2020). Our dynamical models explore different parameterizations of the gravitational potential and are based on two or three distinct stellar populations, separated by metallicity and age.

The article is organized as follows: Section 3 describes the MUSE data set and briefly summarizes the relevant results of the first two papers in this series; Section 4 describes and presents the results of our discrete population-dynamical models; Section 5 discusses the results in the context of previous work; and Section 5 summarizes our findings.

For consistency, we take the base parameters for M54 from the Harris catalog of GCs (Harris 1996, 2010 edition). They are summarized in Table 1.

## 2. Observations

In this section, we present the three data sets from MUSE (Bacon et al. 2014), which we use to extract individual stellar spectra and create the M54 stellar catalog. MUSE is located at UT 4 ( Yepun) of the VLT at the Paranal Observatory in Chile.

### 2.1. MUSE Wide Field Mode

One part of the data set used to analyze the dynamical properties of M54 consists of our MUSE mosaic observations of the system taken between 2015 June 29 and September 19 during run 095.B-0585(A) (P.I.: Lützgendorf). These observations and the data reduction are described in depth in Paper I and Paper II. Here, we briefly lay out the more important details. We have a 4 × 4 pointing MUSE wide field mode (WFM) seeing-limited mosaic, centered on M54. The field of
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Color images obtained from the MUSE cubes using synthetic *i*, *r*, and *z* filters. The left panel shows the WFM MUSE data composed of a 4 × 4 mosaic of 16 pointings, covering an area of ~2.5 $r_h$ (the same as in Paper I; $r_h = 49''$, Harris 1996, 2010 edition). The blue square shows the position of the single WFM+AO 1' × 1' pointing. The WFM+AO color image is presented in the top right panel. The orange square shows the position of the single NFM pointing, which color image is presented in the bottom right panel with a corresponding size of 7''5 × 7''5. North is up and east is to the left.</figcaption></figure>

view of a single MUSE pointing is 59''9 × 60''0 with a
wavelength coverage in the range of 4800–9300 Å, and a
spatial sampling of 0''2 pix⁻¹. Due to some overlap between
the individual pointings, our data covers a square field of view
with a side of 210'' (27.3 pc, assuming a 26.5 kpc distance; see
Table 1). Figure 1 shows a synthetic color image of the MUSE
mosaic, constructed from the MUSE IFU cubes. According to
the GC catalog of Harris (1996, 2010 edition), M54 has a core
radius $r_c = 5''4$ (0.7 pc), a half-light radius $r_h = 49''$ (6.4 pc),
and a tidal radius $r_t = 592''$ (77 pc). Thus, our data covers the
system out to $\gtrsim 2$ half-light radii.

## 2.2. MUSE WFM with Adaptive Optics

We observed adaptive optics (AO) corrected WFM pointing
(WFM+AO), centered on M54 on 2017 August 8 with the
MUSE AO facility, which utilizes the GALACS1 AO module
for ground layer atmospheric turbulence correction, as part of
the MUSE WFM+AO science verification program 60.A-9181
(A) (PI: Alfaro-Cuello). The WFM+AO configuration is the
same as the seeing-limited WFM observations (the same
wavelength coverage and spectral and spatial resolution);
however, due to the sodium lasers used by the AO facility that
saturate the detector, there is a gap in the spectra between 5780
and 6050 Å. The point-spread function (PSF) FWHM of the
WFM+AO field is 0''6. We took four science exposures (750 s
each), applying 90° field rotation between them without
dithering. We show the color image obtained from the WFM

+AO data cube using the synthetic *i*, *r*, and *z* filters in the top right panel of Figure 1.

## 2.3. MUSE Narrow Field Mode with AO

We also include in this analysis a central narrow field mode
(NFM+AO) data cube in the very center of M54, taken during
the MUSE/NFM science verification observations on 2018
September 8 (program ID 60.A-9486(A), P.I. Alfaro-Cuello,
first results shown in Leibundgut et al. 2019). MUSE NFM
+AO mode utilizes a laser tomographic AO system that pro-
vides a 7''42 × 7''43 field of view and nearly diffraction-
limited images. This single central pointing covers the inner-
most region of M54 corresponding to roughly 1 × 1 pc at the
assumed distance of 26.5 kpc. The PSF FWHM is better than
0''2. We present the color image obtained from the NFM cube
using synthetic *i*, *r*, and *z* filters in the bottom right panel of
Figure 1. The MUSE NFM mode has the same wavelength
coverage as the WFM (4800–9300 Å), but with a 10 times
better spatial sampling of 0''025 pix⁻¹. Similar to the WFM
+AO observations, there is a gap in the spectra between 5780
and 6050 Å due to light contamination from the sodium lasers.
We acquired four science exposures (900 s each) and no rota-
tion nor dithering was applied between exposures.

This latter data set is especially important to resolve indi-
vidual stellar spectra in the innermost region of the cluster and
will be crucial to constrain the central mass profile and for a
possible detection of an intermediate-mass black hole in the
center of Sgr's NSC, where crowding severely limits the non-AO data (M. Alfaro-Cuello et al., in preparation).

## 2.4. Individual Stellar Spectra Extraction and Analysis

The stellar spectra of resolved stars were extracted with PAMPELMUSE (Kamann et al. 2013).<fnref n="12" /> This software models the change in the point-spread function with wavelength, allowing to de-blend and extract the spectra of sources efficiently even in crowded and dense regions. This program needs a photometric reference catalog, for which we use Hubble Space Telescope (HST) photometry from Sarajedini et al. (2007), described in Siegel et al. (2007). However, as we pointed out in Paper II, even with PAMPELMUSE we were not able to fully resolve the spectra of the innermost stars, which led to an artificial drop of the velocity dispersion in the central region. This problem has been significantly alleviated by including the MUSE WFM-AO and NFM observations.

During the stellar spectra extraction with PAMPELMUSE we model the PSF using a Moffat profile for the WFM and WFM+AO MUSE data sets. For the NFM data set, we use the MAOPPY PSF model profile (Fétick et al. 2019), which is better suited better for AO observations and especially the MUSE NFM cubes.

Our final stellar catalog in the field of view of M54 contains 8927 entries with measured radial velocities and metallicities from spectra with signal-to-noise ratios better than 10 px⁻¹. Of these, 3008 sources are extracted from the WFM-AO cube and 627 from the NFM cube. We analyzed all spectra from the three data cubes in a uniform manner with the software tool SPEXXY<fnref n="13" />,<fnref n="14" /> (Husser et al. 2016). SPEXXY is a full spectrum fitting framework written in PYTHON with the purpose to derive the main stellar parameters from an observed spectrum by fitting a grid of model spectra. For this work, we utilized the PHOENIX library of synthetic spectra<fnref n="15" /> (Husser et al. 2013), computed for solar-scaled chemical composition and convolved with the MUSE line spread function, kindly provided by Tim-Oliver Husser. With this setup we derived the radial velocities, metallicities, and effective temperatures for all of our good quality sources, while keeping the stellar surface gravities fixed from expectations according to the Dartmouth set of isochrones<fnref n="16" /> (Dotter et al. 2008). The median radial velocity uncertainty for individual stars is 2.8 km s⁻¹ and the median [Fe/H] uncertainty is 0.08 dex. The radial velocity estimates are corrected for perspective effects using the equations from van de Ven et al. (2006), which are within ±0.2 km s⁻¹ for our field of view. We took M54's bulk proper motion from Vasiliev & Baumgardt (2021; PM<sub>R.A.</sub> = −2.683 ± 0.025 mas yr⁻¹, PM<sub>decl.</sub> = −1.385 ± 0.025 mas yr⁻¹).

We note that there is a significant difference in the setup for deriving stellar metallicities used in this work and in Paper I and hence there are inevitable systematic differences. In this work, we use a different fitting software and adopt a different spectral library with significantly larger wavelength coverage. We also mask out the Na doublet at 5892 Å, which was included in the metallicity fits in Paper I. Most notably, on the MP end of the metallicity distribution, we measure consistently

higher metallicities with respect to Paper I in the order of 0.2 dex. Although both methodologies relied on solar-scaled models, the current models extend significantly further to the red and include the region of the Ca triplet, which can be heavily affected by nonsolar [α/Fe] abundances. Thus, the systematically higher metallicities derived here are not surprising, as MP M54 stars are expected to be α-enhanced. At the metal-rich (MR) end of the metallicity distribution, both methods yield consistent results.

The complete data set used in this work (Appendix B, Table 7) is available online.

## 3. Discrete Population-dynamical Models

In the dynamical analysis of M54, we utilize the Jeans Anisotropic Multi-Gaussian Expansion (JAM) method initially developed by Cappellari (2008). The Jeans equations (Jeans 1915) are derived from the steady-state collisionless Boltzmann equation and describe statistically the motion of a large collection of stars in a gravitational potential. Under the assumption of axisymmetry, there are two equations, which link the functional form of the gravitational potential to the stellar density distribution of the system and velocity moments of the stars.

In our work, we use the PYTHON version of the axisymmetric JAM code, written by M. Cappellari<fnref n="17" /> (Cappellari 2008, 2012). According to Paper II, M54's rotation axis, which is also the symmetry axis in our models, is coincidentally oriented along the sky meridian and thus vertically on the MUSE field of view. We also assume that we see M54 edge-on, based on the 3D velocity analysis with Gaia proper motions published in Paper II.

We build multi-population models with discrete posterior distribution functions (Watkins et al. 2013), following an approach similar to Zhu et al. (2016) and Kamann et al. (2020), which describe the distinct kinematic and morphological signatures of multiple stellar populations self-consistently. The models are fully probabilistic and follow a Bayesian framework, i.e., we estimate the probability of each star from our observational sample belonging to a stellar population (k), based on its coordinates, radial velocity, age, and metallicity, given the respective uncertainties of these quantities. There are three ingredients that describe the different stellar populations k and constitute our dynamical models, which we describe as probability density functions (PDF):

1. $P_{\text{pop}}^k$—a stellar population PDF, based on the inferred expected metallicity (and age) distributions of each population. In the two-population models, $P_{\text{pop}}^k$ is Gaussian, characterized by a mean metallicity and a metallicity spread. In the three-population model, $P_{\text{pop}}^k$ is a multivariate Gaussian, characterized by a mean metallicity, a metallicity spread, and mean F606W and F814W magnitudes per star, given an isochrone model, and the respective photometric color spread (see Section 3.3 for more details).

2. $P_{\text{spa}}^k$—the spatial distribution or surface luminosity density of each population, which is inferred from the solution of the Jeans equations under the condition that the total surface luminosity density of all populations ($\sum_k P_{\text{spa}}^k$+ foreground) is equal to the observed one.

<footnote n="12">¹² https://pampelmuse.readthedocs.io</footnote>

<footnote n="13">¹³ https://spexxy.readthedocs.io</footnote>

<footnote n="14">¹⁴ https://github.com/thusser/spexxy</footnote>

<footnote n="15">¹⁵ http://phoenix.astro.physik.uni-goettingen.de</footnote>

<footnote n="16">¹⁶ http://stellar.dartmouth.edu</footnote>

<footnote n="17">¹⁷ https://www-astro.physics.ox.ac.uk/~cappellari/software/#jam</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** MGE fit (solid red line) to the observed surface brightness profile of M54 from (Noyola & Gebhardt 2006, thick black line). The individual Gaussian components are plotted with solid blue lines.</figcaption></figure>

3. $P_{\text{dyn}}^k$—each population's velocity moments distribution—also inferred from the solution of the Jeans equations. Assuming a certain shape of the gravitational potential and the population's angular momentum ($\kappa^k$) and anisotropy ($\beta_z^k$), the JAM code predicts the first ($\bar{v}$) and second ($\bar{v}^2$) velocity moments at the position of each star in the sample. We compute the probability of each star's radial velocity to be drawn from a Gaussian with the first line-of-sight velocity moment as the mean and a variance $\sigma^2 = \bar{v}^2 - \bar{v}^2$. The anisotropy is defined as $\beta_z = 1 - \frac{\bar{v}^2}{\bar{v}^2 - \bar{v}^2}$.

where $\frac{\bar{v}^2}{\bar{v}^2 - \bar{v}^2}$ is the ratio of the second velocity moments in the vertical and radial direction in cylindrical coordinates. The angular momentum parameter $\kappa$ is a dimensionless quantity that provides a direct measure of the population's rotation, defined as $\kappa = \frac{\bar{v}_\phi}{(\bar{v}_\phi^2 - \bar{v}^2)^{1/2}}$, where $\bar{v}_\phi$ and $\bar{v}^2$ are the first and second velocity moments in the tangential direction. When $\kappa = \pm 1$ and $\beta_z = 0$, the system reduces to an isotropic rotator and when $\kappa = 0$, there is no net angular momentum. Negative $\beta_z$ values indicate tangential anisotropy, while positive $\beta_z$ values indicate radial anisotropy (Cappellari 2008). Although both the anisotropy and rotation parameters can vary with radius in GCs (Bianchini et al. 2013, 2017a), in this work we assume that they are radially constant, but independent for each population.

The JAM code uses multi-Gaussian expansion (MGE) models to describe the surface luminosity density and the gravitational potential (Emsellem et al. 1994). We use the V-band surface brightness profile of Noyola & Gebhardt (2006) as representative of the system and fit it with a 10-component MGE model (Figure 2), using the MGE implementation and software of Cappellari (2002). According to Harris' catalog (2010 version) the integral magnitude of M54 is $M_V = -9.98$ mag, which corresponds to an integral luminosity of $L_V = 0.85 \times 10^6 L_\odot$. We scale the MGE to match this figure.

We also include in all models an additional population component with a flat surface density, which represents the foreground stellar distribution. Its radial velocity moments are derived from the Besançon model of the Galaxy (Robin et al. 2003) at the position of M54, which is integrated along the line

of sight. ($\langle V \rangle = 40 \text{ km s}^{-1}$, $\sigma_V = 90 \text{ km s}^{-1}$). The foreground population is expressed with a single free parameter ($\epsilon$), as a fraction of the central surface brightness of M54, as in Watkins et al. (2013) and Zhu et al. (2016). Hence, we can write the spatial PDFs of the different model components as

$$P_{\text{spa}}^k = S_k / \sum_k (S_k + \epsilon C_k) \quad (1)$$

and

$$P_{\text{spa}}^{\text{foreground}} = 1 - \sum_{k} P_{\text{spa}}^{k}, \quad (2)$$

where $S_k$ is the surface brightness profile and $C_k$—the central surface brightness of the various modeled populations ($k$).

The posterior probability of a star ($i$) belonging to either of the dynamical components or to the foreground component of the model ($k$) is then given by the joint probability of the above specified three probability distributions and we maximize the log posterior function of the entire sample:

$$\ln P = \sum_{i} \ln \sum_{k} P_{\text{pop},i}^{k} P_{\text{spa},i}^{k} P_{\text{dyn},i}^{k} + \ln P_{\text{pri}}, \quad (3)$$

where $P_{\text{pri}}$ contains the parameter priors (see below). We use the EMCEE affine invariant Markov Chain Monte Carlo (MCMC) algorithm (Goodman & Weare 2010) in PYTHON (Foreman-Mackey et al. 2013). We run our models on a CPU cluster engaging 96 cores and use 192 walkers and 6000 steps in the MCMC, which we confirmed to be enough for the fit to converge.

We explore dynamical models with two- and three-population components. In the two-population model, we use only the measured metallicity of each star in the sample as a population tag, while in the three-population model we also include the F606W and F814W HST magnitudes and fit for the population ages using isochrones. The former model is simpler in construction and it aims to describe the interplay between the MP and MR stars in the M54 system, but ignores the Sgr's field population, which is largely mixed with the MR nuclear component, but has different spatial and kinematic properties. The three-population model, on the other hand, aims to provide a more complete physical picture of this dynamically complex system, exploring the interplay of several population-dynamical components simultaneously, but requires more prior assumptions to converge to a physical solution.

### 3.1. Two-population Model

At first, we consider a population-dynamical model with two distinct stellar populations, separated by metallicity. We assume that M54 consists of an MP and MR component and we fit for their mean metallicities and intrinsic metallicity spreads, which gives four free parameters.

To model the spatial probability distributions of the two populations, we allow them to have decoupled surface brightness profiles, by setting each population to contribute an independent fraction to the luminosity of each of the 10 MGE components, which define the surface brightness of the entire nuclear system. This approach ensures that the sum of the two decoupled surface brightness profiles is equal to the observed one from Noyola & Gebhardt (2006). In principle, the JAM models can also take a flattening value ($q = \frac{b}{a}$, where $a$ and $b$
are the major and minor axes, respectively) for each MGE component, but here we adopt a radially constant flattening for each population, and assume that the two flattening parameters are independent. This leads to 12 additional free model parameters—one fraction for each MGE component and two flattening values: $q_{MP}$ and $q_{MR}$ for the MP and MR populations, respectively.

For the gravitational potential we first assume that mass follows light and introduce a constant mass-to-light ratio ($M/L_V \equiv \Upsilon_V$) for each of the two populations as free parameters. We impose weak Gaussian priors with relatively large variance on these quantities ($\langle\Upsilon_V^{\text{MP}}\rangle = 2.2 M_\odot/L_\odot$, $\sigma_\Upsilon^{\text{MP}} = 0.5 M_\odot/L_\odot$ and $\langle\Upsilon_V^{\text{MR}}\rangle = 1.1 M_\odot/L_\odot$, $\sigma_\Upsilon^{\text{MR}} = 0.5 M_\odot/L_\odot$) so that they are loosely consistent with expectations from stellar population synthesis models. According to the BaSTI<fnref n="18" /> stellar population models (Percival et al. 2009), an old MP population (age = 11.5 Gyr, [Fe/H] = -1.27) is expected to have $\Upsilon_V = 2.2 M_\odot/L_\odot$ and a YMR population (age = 2.5 Gyr, [Fe/H] = -0.25) is expected to have $\Upsilon_V = 1.1 M_\odot/L_\odot$. The total mass profile of M54 can then be computed in the following way:

$$M(r) = \Upsilon_{\text{MR}} L_{\text{MR}}(r) + \Upsilon_{\text{MP}} L_{\text{MP}}(r). \quad (4)$$

We explore a gravitational potential with radially variable $\Upsilon(r)$ later in Section 3.2.

In addition, we set as free parameters dimensionless quantities that describe each population's angular momentum ($\kappa_{MR}$, $\kappa_{MP}$) and velocity anisotropy ($\beta_{z,MR}$, $\beta_{z,MP}$). We also set the systemic velocity (RV) of the system as a free parameter in this model ($\langle V \rangle$). An additional parameter controls the foreground fraction of stars ($\epsilon$). This leaves us with a total of 24 free model parameters.

A corner plot showing the posterior distributions of the 14 main parameters is presented in Figure 3. Their values are summarized in Table 2. In Appendix A (Figure 17) we also show the covariances of the 10 additional parameters describing the MGE fractions contributing to the MP population (the MR population fractions are one minus those, respectively). One could note from Figure 17 that not all MGE components are relevant for separating the MP and MR populations' surface brightness distributions; however, their combined effect, marginalized over the posterior, leads to a reliable and reproducible solution.

We find dynamical $\Upsilon_V^{\text{MP}} = 2.3 \pm 0.1 M_\odot/L_{\odot,V}$ and $\Upsilon_V^{\text{MR}} = 0.2 \pm 0.2 M_\odot/L_{\odot,V}$. $\Upsilon_V^{\text{MP}}$ is consistent with expectations from synthesis models of old (12 Gyr) MP stellar populations at this metallicity according to the BaSTI stellar population synthesis tool. $\Upsilon_V$ for the MR stars of M54 on the other hand is significantly lower than expectations for 1.5–2.0 Gyr MR ([Fe/H] = -0.25 dex) populations synthesis models and consistent with zero. The model tends to attribute all mass to only one of the two populations.

With this perhaps overly simplistic assumption of constant $\Upsilon$ ratios for the two stellar populations, our model predicts a total mass of the system of $1.26 \pm 0.03 \times 10^6 M_\odot$ at the tidal radius. The mass of the MP population is 96% of the total, while it accounts for only 64% of the luminosity. As a result, we estimate a global $\Upsilon_V = 1.51 \pm 0.04 M_\odot/L_\odot$ for the entire system, which is in line with expectations for a mixed system consisting of both old and young stellar populations.

We also ran a model with a single constant $M/L$ ratio describing the gravitational potential for the entire system. Not surprisingly, we find a $\Upsilon_V = 1.52 \pm 0.03 M_\odot/L_{\odot,V}$, which is very similar to the global value found above. We note that in this test we did not impose any prior on $\Upsilon$.

When we refer to the total dynamical mass in this work, we consider the mass locked within the tidal radius of the system, as given in the Harris catalog of GCs (see Table 1; Harris 1996, 2010 edition), which corresponds to a radial distance of ~77 pc. However, these values need to be taken cautiously for two main reasons. First, our kinematic data does not reach that far out in the system and therefore, the quoted quantities are extrapolated predictions of the model. Our MUSE mosaic field of view contains ~85% of the total luminosity of the nuclear system. Second, despite feeling the Milky Way tidal field, M54's nucleus is still embedded in the densest central region of the disrupting Sgr dSph and hence its tidal radius is not well defined. The tidal radius that we quote comes simply from the best-fit King (1962) profile to the surface brightness profile of the system by McLaughlin & van der Marel (2005). We stress that the surface brightness profile that we use in our dynamical models (approximated with an MGE) is not truncated at this radius. In this respect, it is much more meaningful to look at the cumulative mass profile that we obtain with our dynamical models and present in Figure 4. We present the estimated masses of the two populations and the entire system at several radial distances in Table 3.

We conclude that with this choice of gravitational potential parameterization, the population-dynamical model is prone to attributing all of the mass to a single component and hence we next explore dynamical models with a single radially varying $\Upsilon$. To avoid repetition, we discuss the other results of the two-population-dynamical model and their implications in the next section, as they appear practically unaffected by the choice of the gravitational potential parameterization.

## 3.2. Radially Varying M/L Ratio

The population-dynamical model discussed in the previous section has one significant drawback, which is the tendency to attribute the entire mass budget to only one of the model components. Our solution is to model the gravitational potential globally for the entire system, but allow for more degrees of freedom when fitting for it, by introducing a radially varying universal $M/L$ ratio, $\Upsilon(r)$.

We adopted a technique similar to that of Kamann et al. (2020) to model the gravitational potential. Essentially, we assigned an individual $\Upsilon$ to each MGE component in the surface brightness profile of M54, according to a parameterization given by Equation (5) in Kamann et al. (2020):

$$\Upsilon(r) = \frac{\Upsilon_0 \left(1 - \frac{r}{r_{\Upsilon}}\right) + 2\Upsilon_1 \left(\frac{r}{r_{\Upsilon}}\right) - \Upsilon_{\infty} \frac{r}{r_{\Upsilon}} \left(1 - \frac{r}{r_{\Upsilon}}\right)}{1 + \left(\frac{r}{r_{\Upsilon}}\right)^2}, \quad (5)$$

where $\Upsilon_0$ is the $M/L$ ratio in the center of the system, $\Upsilon_\infty$, at infinity, and $\Upsilon_t$ is a transition value at a characteristic radius, $r_\Upsilon$. The $\Upsilon$ value corresponding to each MGE component is calculated at a radius of $1\sigma$ of that component, where $\sigma$ is the standard deviation width of each Gaussian component, comprising the MGE. Hence, the resulting $M/L_V$ profile from dividing the mass $M(r)$ and luminosity $L_V(r)$ profiles from the MGE components is

<footnote n="18">¹⁸ http://albione.oa-teramo.inaf.it</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** MCMC corner plot showing the posterior distributions of the two-population-dynamical model parameters with two constant *M/L* ratios (Section 3.1). The 10 light fractions of the surface brightness MGE are omitted.</figcaption></figure>

not exactly the same as $\Upsilon(r)$ in Equation (5). The functional form of the radially varying $\Upsilon$ is motivated by GC dynamical models (den Brok et al. 2014a; Bianchini et al. 2017b) and N-body simulations (Baumgardt 2017), which show that the *M/L* ratio radial profile of GCs is cup-shaped with a peak in the center, a gradual decrease toward the half-mass-radius, followed by another increase in the outer regions. Equation (5) mimics this behavior. These three $\Upsilon$ figures and the characteristic radius are free parameters in our model. Baumgardt (2017) note that for the majority of Galactic GCs the characteristic radius ($r_\Upsilon$) is at about a tenth of the half-mass-radius ($r_h \simeq 50''$ or 6.8 pc), so we introduce a Gaussian prior with this condition ($\langle r_\Upsilon \rangle = 5''$, $\sigma_{r_\Upsilon} = 1''$), otherwise $\Upsilon_t$ and $r_\Upsilon$ are degenerate. Furthermore, as our data is spatially limited and does not extend to the tidal

radius of M54, we do not have a good constraint on $\Upsilon_\infty$, so we also set a Gaussian prior on this parameter with a mean of $\Upsilon_\infty = 2 M_\odot/L_\odot$ and a standard deviation 0.5 $M_\odot/L_\odot$, which is typical for GCs and ultracompact dwarfs old stellar populations (Voggel et al. 2019).

We kept all other model parameters the same as in Section 3.1, which results in a population-dynamical model with 26 free model parameters. A corner plot showing their covariances is presented in Figure 5 and their values are summarized in Table 2.

Figure 6 shows the resulting radial $\Upsilon_V$ ratio profile and the estimated cumulative mass distribution of M54. The best-fit model indeed prefers a slightly higher $\Upsilon$ in the innermost region of the system with a minimum at the core radius, followed by a gradual increase outwards. It is assumed that the
<caption>Table 2</caption>

Dynamical Model Best-fit Parameters—Median and Standard Deviations from the MCMC Posterior Distribution

<table><thead><tr><th></th><th>Section 3.1</th><th>Section 3.2</th><th>Section 3.3</th></tr></thead><tbody><tr><td>⟨V<sub>0</sub>⟩ (km s<sup>-1</sup>)</td><td>142.50 ± 0.13</td><td>142.46 ± 0.13</td><td>...</td></tr><tr><td>Υ<sub>MR</sub>(M<sub>☉</sub>/L<sub>☉</sub>,V)</td><td>0.20 ± 0.18</td><td>...</td><td>...</td></tr><tr><td>Υ<sub>MP</sub>(M<sub>☉</sub>/L<sub>☉</sub>,V)</td><td>2.27 ± 0.11</td><td>...</td><td>...</td></tr><tr><td>Υ<sub>0</sub>(M<sub>☉</sub>/L<sub>☉</sub>,V)</td><td>...</td><td>5.42 ± 2.62</td><td>5.62 ± 2.35</td></tr><tr><td>Υ<sub>f</sub>(M<sub>☉</sub>/L<sub>☉</sub>,V)</td><td>...</td><td>1.10 ± 0.25</td><td>1.00 ± 0.15</td></tr><tr><td>Υ<sub>∞</sub>(M<sub>☉</sub>/L<sub>☉</sub>,V)</td><td>...</td><td>2.89 ± 0.29</td><td>3.15 ± 0.36</td></tr><tr><td>r<sub>f</sub>(")</td><td>...</td><td>3.98 ± 1.25</td><td>4.38 ± 0.80</td></tr><tr><td>q<sub>YMR</sub></td><td>0.90 ± 0.04</td><td>0.91 ± 0.03</td><td>0.84 ± 0.03</td></tr><tr><td>q<sub>IMR</sub></td><td>...</td><td>...</td><td>0.88 ± 0.06</td></tr><tr><td>q<sub>OMP</sub></td><td>0.96 ± 0.03</td><td>0.97 ± 0.03</td><td>0.95 ± 0.04</td></tr><tr><td>β<sup>YMR</sup></td><td>-0.09 ± 0.06</td><td>-0.08 ± 0.06</td><td>0.00 ± 0.08</td></tr><tr><td>β<sup>IMR</sup></td><td>...</td><td>...</td><td>0.09 ± 0.14</td></tr><tr><td>β<sup>OMP</sup></td><td>0.04 ± 0.04</td><td>0.04 ± 0.03</td><td>0.08 ± 0.04</td></tr><tr><td>κ<sub>YMR</sub></td><td>-0.63 ± 0.08</td><td>-0.65 ± 0.08</td><td>-0.98 ± 0.17</td></tr><tr><td>κ<sub>IMR</sub></td><td>...</td><td>...</td><td>-0.06 ± 0.34</td></tr><tr><td>κ<sub>OMP</sub></td><td>-0.52 ± 0.53</td><td>-0.31 ± 0.52</td><td>0.07 ± 1.05</td></tr><tr><td>⟨[Fe/H]⟩<sub>YMR</sub></td><td>-0.28 ± 0.007</td><td>-0.28 ± 0.007</td><td>-0.21 ± 0.004</td></tr><tr><td>⟨[Fe/H]⟩<sub>IMR</sub></td><td>...</td><td>...</td><td>-0.47 ± 0.006</td></tr><tr><td>⟨[Fe/H]⟩<sub>OMP</sub></td><td>-1.21 ± 0.003</td><td>-1.21 ± 0.003</td><td>-1.21 ± 0.004</td></tr><tr><td>σ<sup>YMR</sup><sub>[Fe/H]</sub></td><td>0.21 ± 0.006</td><td>0.21 ± 0.006</td><td>0.17 ± 0.008</td></tr><tr><td>σ<sup>IMR</sup><sub>[Fe/H]</sub></td><td>...</td><td>...</td><td>0.26 ± 0.02</td></tr><tr><td>σ<sup>OMP</sup><sub>[Fe/H]</sub></td><td>0.21 ± 0.003</td><td>0.21 ± 0.003</td><td>0.24 ± 0.005</td></tr><tr><td>⟨age⟩<sub>YMR</sub>(Gyr)</td><td>...</td><td>...</td><td>2.53 ± 0.003</td></tr><tr><td>⟨age⟩<sub>IMR</sub>(Gyr)</td><td>...</td><td>...</td><td>5.37 ± 0.06</td></tr><tr><td>⟨age⟩<sub>OMP</sub>(Gyr)</td><td>...</td><td>...</td><td>14.14 ± 0.13</td></tr><tr><td>σ<sup>YMR</sup><sub>V-I</sub>(mag)</td><td>...</td><td>...</td><td>0.05 ± 0.01</td></tr><tr><td>σ<sup>IMR</sup><sub>V-I</sub>(mag)</td><td>...</td><td>...</td><td>0.0005 ± 0.0004</td></tr><tr><td>σ<sup>OMP</sup><sub>V-I</sub>(mag)</td><td>...</td><td>...</td><td>0.0003 ± 0.0003</td></tr><tr><td>E(B-V) (mag)</td><td>...</td><td>...</td><td>0.11 ± 0.0005</td></tr><tr><td>ε [×10<sup>-4</sup>]</td><td>2.6 ± 0.15</td><td>2.7 ± 0.15</td><td>3.1 ± 0.15</td></tr></tbody></table>

mass segregation of stellar remnants is the reason for the upturn in M/L at the center (Bianchini et al. 2017b), while Υ is high in the outer parts due to the overdensity of lower mass stars. In the more complicated case of M54, the presence of the more centrally concentrated MR, significantly younger population (expected lower Υ<sub>V</sub>) works to mitigate the effect of mass segregation. At the same time, the entire system is likely embedded in a dark matter halo (Carlberg & Grillmair 2022), not surprising for a dwarf galaxy nucleus, which further increases the M/L ratio in the outer regions.

Due to the change in the shape of the potential, compared to the dynamical model described in Section 3.1, we find a higher mass for the M54 system (M<sub>tot</sub> = 1.58 ± 0.07 × 10<sup>6</sup> M<sub>☉</sub> at the tidal radius), where 67% of the mass and luminosity are attributed to the MP population. We note that since in this dynamical model we use a common M/L ratio for both populations, the mass contribution of each population mostly follows its luminosity contribution. Small differences arise due to the different spatial extent of the luminosity density distributions of the two populations. We find a global Υ<sub>V</sub> = 1.92 ± 0.08 M<sub>☉</sub>/L<sub>☉</sub> for the entire system, using a two-population-dynamical model with a radially variable M/L ratio. Table 4 gives the mass estimates and Υ<sub>V</sub> of the two populations locked within several key radial points.

We can also immediately see from Figures 3 and 5 that the MP population is close to isotropic (β<sub>z</sub> = 0.04 ± 0.04) with considerably lower angular momentum than the MR one. The

<figure><img src="image_7.png" /><figcaption>Figure 4. Cumulative masses of the two populations (blue—MP; red—MR) and the entire system (black), according to the model described in Section 3.1. The three black vertical dashed lines on both panels correspond to M54’s core, half-light, and tidal radii (GC catalog of Harris (1996), 2010 edition), while the red dashed line indicates the end of the MUSE mosaic field of view.</figcaption></figure>

latter also has a negative β<sup>z</sup> = -0.09 ± 0.06 parameter, which is indicative of a slight tangential anisotropy (Cappellari 2008).

The discrete Jeans model predicts the velocity moments at the position of each observed star and the parameter optimization was performed on the discrete data. To visually compare the model prediction to the observed kinematics, we assign a parent population to each of the observed stars based on a 50% probability threshold, which is only relevant for presentation purposes and plotting. Then we bin both the observations and the models in the same way to produce radial rotation and velocity dispersion profiles (Figure 7). The choice of binning is purely for visualization purposes and does not have any effect on the results and conclusions in this work. We obtain the mean velocity and intrinsic velocity dispersion for each bin in Figure 7 via a maximum likelihood fit to the kinematic distribution using a Gaussian model with unknown mean and variance. The variance in the likelihood function includes the individual measurement uncertainties. We also correct the observed velocity dispersion profile for the systemic rotation. For the rotation profile, we use overlapping linearly spaced vertical slices along the horizontal x-axis (Figure 7, left panel) and for the velocity dispersion, we use logarithmically spaced concentric radial bins (Figure 7, right panel). We perform multiple random draws from the model posterior to derive the model uncertainties. Overall, the two-component Jeans dynamical model describes the observed kinematics of the system at all radii very well. It is evident that the MR population has a high degree of intrinsic rotation combined with a lower velocity dispersion, while the MP component has a significantly lower angular momentum and higher velocity dispersion.

The radial distributions of the two populations, together with their cumulative distribution functions are presented in Figure 8. They are both very close to spherical with axis ratios q<sub>MR</sub> = 0.90 ± 0.04 and q<sub>MP</sub> = 0.96 ± 0.03 for the MR and MP populations, respectively. Note that we reported larger flattening for both populations in Paper I, using a direct elliptical Plummer model fit to their observed stellar number densities - q<sub>MR</sub> = 0.69 ± 0.10 and q<sub>MP</sub> = 0.84 ± 0.06. The method used in Paper I, however, is sensitive to incompleteness in the data, while here we infer the population flattenings purely from
<figure><img src="image_8.png" /><figcaption>**Figure 5.** MCMC corner plot showing the posterior distributions of the two-population-dynamical model parameters with a radially varying M/L ratio (Section 3.2). The 10 light fractions of the surface brightness MGE are omitted.</figcaption></figure>

<table><caption>**Table 3**
Cumulative Dynamical Mass Estimates and $\Upsilon_V$ of M54's Populations According to the Model Presented in Section 3.1</caption><thead><tr><th></th><th>MP<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>MR<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>Total<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>% MP</th><th>% MR</th><th>Υ<sub>MP</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th><th>Υ<sub>MR</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th><th>Υ<sub>Tot</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th></tr></thead><tbody><tr><td>r<sub>c</sub></td><td>0.10 ± 0.01</td><td>0.01 ± 0.01</td><td>0.11 ± 0.01</td><td>92</td><td>8</td><td>2.27 ± 0.08</td><td>0.20 ± 0.17</td><td>1.20 ± 0.05</td></tr><tr><td>r<sub>h</sub></td><td>0.82 ± 0.02</td><td>0.03 ± 0.03</td><td>0.85 ± 0.02</td><td>96</td><td>4</td><td>2.27 ± 0.04</td><td>0.20 ± 0.19</td><td>1.53 ± 0.04</td></tr><tr><td>r<sub>fov</sub></td><td>1.10 ± 0.05</td><td>0.04 ± 0.04</td><td>1.14 ± 0.03</td><td>96</td><td>4</td><td>2.27 ± 0.04</td><td>0.20 ± 0.19</td><td>1.55 ± 0.04</td></tr><tr><td>r<sub>t</sub></td><td>1.21 ± 0.06</td><td>0.05 ± 0.05</td><td>1.26 ± 0.03</td><td>96</td><td>4</td><td>2.27 ± 0.06</td><td>0.20 ± 0.19</td><td>1.51 ± 0.04</td></tr></tbody></table>

Note. The various radii are defined in Table 1.
<table><caption>**Table 4**
Cumulative Dynamical Mass Estimates and $\Upsilon_V$ of M54's Populations According to the Model Presented in Section 3.2</caption><thead><tr><th></th><th>MP<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>MR<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>Total<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>% MP</th><th>% MR</th><th>Υ<sub>MP</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th><th>Υ<sub>MR</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th><th>Υ<sub>Tot</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th></tr></thead><tbody><tr><td>r<sub>c</sub></td><td>0.06 ± 0.01</td><td>0.06 ± 0.01</td><td>0.12 ± 0.01</td><td>50</td><td>50</td><td>1.32 ± 0.07</td><td>1.36 ± 0.12</td><td>1.35 ± 0.14</td></tr><tr><td>r<sub>h</sub></td><td>0.61 ± 0.02</td><td>0.26 ± 0.01</td><td>0.87 ± 0.02</td><td>70</td><td>30</td><td>1.65 ± 0.03</td><td>1.51 ± 0.06</td><td>1.59 ± 0.04</td></tr><tr><td>r<sub>fov</sub></td><td>0.93 ± 0.03</td><td>0.38 ± 0.02</td><td>1.31 ± 0.04</td><td>71</td><td>29</td><td>1.85 ± 0.03</td><td>1.71 ± 0.06</td><td>1.81 ± 0.06</td></tr><tr><td>r<sub>t</sub></td><td>1.06 ± 0.04</td><td>0.52 ± 0.04</td><td>1.58 ± 0.07</td><td>67</td><td>33</td><td>1.93 ± 0.05</td><td>1.90 ± 0.11</td><td>1.92 ± 0.08</td></tr></tbody></table>

Note. The various radii are defined in Table 1.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Dynamical results according to the model described in Section 3.2. Top: M54 radial $\Upsilon_V$ profile, resulting from dividing the best-fit mass and luminosity MGE models of the system. The median $\Upsilon_V$ profile is indicated with a thick curve, while the thin curves represent random draws from the posterior distribution. The horizontal dashed line corresponds to the best-fit constant $\Upsilon_V$. Bottom: Cumulative masses of the two populations (blue—MP; red—MR) and the entire system (black). The three black vertical dashed lines on both panels correspond to M54's core, half-light, and tidal radii (GC catalog of Harris (1996, 2010 edition), while the red dashed line indicates the end of the MUSE mosaic field of view.</figcaption></figure>

dynamical constraints. It is important to note here that the modeled density profiles of the two populations are inferred solely from the Jeans equations and the observed stellar kinematics, and do not depend on the observed stars' spatial distribution, which is heavily biased due to the very uneven photometric depth reached by the different MUSE data sets. In fact, there is no need that the observed kinematic tracers follow

the actual density distribution of the stars in order to fit for their density profiles. This makes our models insensitive to various sources of incompleteness in the observations and we can also predict the stellar density outside of the MUSE mosaic field of view, relying on the M54's surface brightness profile by Noyola & Gebhardt (2006), which extends much further out. Figure 8 shows that the MP and MR components have different spatial distributions. MR stars are generally more centrally concentrated than MP stars, but their distribution flattens out significantly in the outer region. We noted that MR stars are already more centrally concentrated than MP stars in Paper I, but here we show that this is a prediction of the best-fit dynamical model too. The excess of MR stars in the outer region belong to intermediate-age MR Sgr field stars, which we do not explicitly consider in this model and they appear mixed with the MR M54 population. In Section 3.3, we introduce a three-component dynamical model to alleviate this problem.

In Figure 9, we show the separation between the two dynamically distinct components using two-population diagnostic plots. The left panel of Figure 9 shows M54's observed stellar sample in radial velocity versus metallicity space, where individual stars are color coded with varying shades of blue and red for the MP and MR populations, respectively, depending on their membership probability. The two populations with identical mean systemic velocities are well distinguishable on the [Fe/H] axis. We find a mean $\langle[Fe/H]_{MR}\rangle = -0.28 \pm 0.01$ dex and $\langle[Fe/H]_{MP}\rangle = -1.21 \pm 0.01$ dex with both populations having the same intrinsic metallicity spread of $0.21 \pm 0.01$ dex. The reported intrinsic metallicity spreads take into account the [Fe/H] measurement uncertainties estimated with the SPEXXY code; however, there are indications that these uncertainties could be underestimated by 20%–30% (Husser et al. 2016). Still the measured intrinsic spreads are quite large and there are likely physical reasons for this. The large [Fe/H] spread in the MR population is likely due to an inhomogeneous origin of the MR stars, which we explore in Section 3.3, while the large [Fe/H] spread in the MP population could be explained by merging GCs with a spread in metallicities.

In the right panel of Figure 9, we show the HST color-magnitude diagram (CMD) of our sample, using photometry from Sarajedini et al. (2007) and Siegel et al. (2007) and the same population probability color coding. As expected, the MP stars follow a narrow distribution on the CMD that is consistent with an old stellar population, similar to the majority of halo GCs. On the other hand, the MR stars are clearly younger and likely have a larger age spread.

The differences in kinematics, spatial distribution, and population properties are all well captured in a single self-consistent dynamical model, based on the Jeans equations. The high angular momentum, lower dispersion, central concentration, and high metallicity of the MR stars all point to an in situ
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Predictions for the rotation (left panel) and the velocity dispersion (right panel) for the two populations drawn from the model posterior (MR—red lines, MP—blue lines) according to the dynamical model described in Section 3.2, compared to the observed binned profiles (red squares and blue dots with error bars, respectively), using stars with >50% probability of belonging to either of the two populations. The horizontal error bars indicate the size of each bin. The vertical dashed line in the right panel shows the border between the NFM and WFM MUSE observations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left panel: radial distribution functions of the two-population-dynamical components of M54 as drawn from the Section 3.2 model posterior (MR—red lines and MP—blue lines). The solid black line shows the combined surface brightness profile. Right panel: radial cumulative distribution function of the two-population-dynamical components of M54.</figcaption></figure>

origin of this structure, where it originated in a rapidly rotating gaseous disk in the center of the system. The MR population appears marginally flatter than the MP population, although both are very close to spherical. There are various mechanisms of disk heating, that could puff the MR stars to an almost spherical distribution. For instance, the disk relaxes and redistributes its angular momentum to the stars in the MP population, becoming more spherical with time (see Paper II; Mastrobuono-Battisti & Perets 2013, 2016). On the contrary, the MP population properties are fully consistent with GC origin, where one or more GCs arrived at the center of Sgr through dynamical friction.

### 3.3. Three-population Model

In Paper I we identified the IMR population as representative of the Sgr field stellar content and showed that it has a more extended SFH and surface density distribution than the YMR and OMP populations. It is also distinct with its flat velocity dispersion profile (Paper II).

Here, we present a three-population-dynamical model of M54, which includes simultaneously all three components identified in Papers I and II—the YMR, IMR, and OMP. In

order to better separate the different populations, we also consider their age difference and fit for the three populations’ mean ages. To this aim, we use the F606W and F814W magnitudes and their uncertainties from the Sarajedini et al. (2007) HST catalog in addition to the measured stellar metallicities from the MUSE spectra and fit three isochrones simultaneously with the velocity moments from the Jeans equations to the three modeled populations.

We opted to work with the scaled-solar Dartmouth set of isochrones and we developed an isochrone interpolator to this aim, that generates an isochrone for any custom combination of age and metallicity in the abovementioned photometric bands. We tested isochrone grids with varying [α/Fe], but found out that the scaled-solar isochrones were a better match to the observed CMD, given the estimated stellar metallicities and fixed distance. We note that while it is expected that the OMP stars are more α-enhanced than the YMR and IMR stars, we measured the metallicities assuming a scaled-Solar spectral library. In addition, we model the horizontal branch (HB) of the OMP population and the red clump (RC) of the YMR population using the Dartmouth Stellar Evolution Database HB/AGB Track Grids. We chose the HB and RC models that best
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Metallicity vs. radial velocity plot of all stars in the sample (left panel) and a CMD (right panel), color coded by the probability of the star belonging to the MR population (shades of red) or the MP population (shades of blue), according to the dynamical model described in Section 3.2. The photometry is from Sarajedini et al. (2007).</figcaption></figure>

represent M54’s CMD and kept them fixed, e.g., independent of the age and metallicity of the fitted isochrones. We kept the assumed distance to M54, which is relevant for both the isochrone and velocity moments fits, fixed to 26.5 kpc, but set the line-of-sight reddening as a free parameter. Overall, we fit for the mean metallicity and mean age of the YMR, IMR, and OMP populations, their corresponding intrinsic metallicity spreads, reddening, and we allow for a larger variance in the F606W–F814W color in addition to the photometric uncertainties. The latter could eventually be interpreted as coming from an intrinsic age spread of the individual populations, but such a connection is not straightforward and we refrain from making it. We assume a multidimensional Gaussian population likelihood for each observed star to belong to either of the three populations, based on its F606W, F814W magnitudes and metallicity, given an interpolated isochrone model defined by age, metallicity, and reddening.

We use the gravitational potential with a radially varying $\Upsilon$, as described in Section 3.2 and adopt the same priors in the posterior function: $\langle r_{\Upsilon} \rangle = 5''$, $\sigma_{r_{\Upsilon}} = 1''$ and $\Upsilon_{\infty} = 2 M_{\odot}/L_{\odot}$ with a standard deviation $0.5 M_{\odot}/L_{\odot}$.

We take a slightly different approach in describing the spatial distribution of the three populations. Instead of letting the Gaussian components in the surface brightness MGE each contribute an independent fraction to each of the three populations (as in Section 3.1), we adopt Legendre polynomials to describe the fractions of each stellar population as a function of radius. While the former approach is maximally agnostic about the shape of the individual populations surface brightness distributions, it requires too many free parameters, not all of which contribute meaningfully to the fit. In a two-population model, we require 10 free parameters, one for each MGE component, in a three-population model, we would require 20 additional free parameters—two for each MGE component, which quickly becomes unfeasible. By choosing Legendre polynomials we greatly reduce the number of additional parameters, at the expense of partly losing the independence of the

individual MGE fractions. We believe, however, that this choice of parameterization is still agnostic enough about the shape of the populations’ surface brightness distributions, so that we do not introduce involuntary bias to the fit. We need two polynomials, as the fraction of the third population is determined by the remaining fraction from the total. Similar to the functional description of the gravitational potential, we use the values of the Legendre polynomials at radii of $1\sigma$ of the MGE components to scale them accordingly. We tested 2nd and 3rd order polynomials and found very similar results. We chose to work with the 2nd order Legendre polynomials for simplicity.

In this model, we also fit for the flattening, angular momentum, and anisotropy of the three populations independently, but we fixed the systemic velocity of M54 to the best-fit value in the two-population model to reduce the number of free parameters.

The three-population model has 33 free parameters: four parameters describing the gravitational potential ($\Upsilon_0$, $\Upsilon_t$, $\Upsilon_\infty$, $r_t$); three anisotropy ($\beta_z^{\text{YMR}}$, $\beta_z^{\text{IMR}}$, $\beta_z^{\text{OMP}}$), three rotation ($\kappa_{\text{YMR}}$, $\kappa_{\text{IMR}}$, $\kappa_{\text{OMP}}$), and three flattening ($q_{\text{YMR}}$, $q_{\text{IMR}}$, $q_{\text{OMP}}$) parameters describing the kinematics of the three populations; three mean ages ($\langle\langle\text{age}\rangle\rangle_{\text{YMR}}$, $\langle\langle\text{age}\rangle\rangle_{\text{IMR}}$, $\langle\langle\text{age}\rangle\rangle_{\text{OMP}}$); three mean metallicities ($\langle\langle[\text{Fe/H}]\rangle\rangle_{\text{YMR}}$, $\langle\langle[\text{Fe/H}]\rangle\rangle_{\text{IMR}}$, $\langle\langle[\text{Fe/H}]\rangle\rangle_{\text{OMP}}$); three intrinsic metallicity spreads ($\sigma_{[\text{Fe/H}]_{\text{IMR}}}^{\text{YMR}}$, $\sigma_{[\text{Fe/H}]_{\text{IMR}}}^{\text{IMR}}$, $\sigma_{[\text{Fe/H}]_{\text{IMR}}}^{\text{OMP}}$); three intrinsic photometric spreads (additional photometric uncertainty $\sigma_{V-I}^{\text{IMP}}$, $\sigma_{V-I}^{\text{IMR}}$, $\sigma_{V-I}^{\text{YMR}}$); one foreground parameter ($\epsilon$); one reddening parameter ($E(B-V)$); three Legendre polynomial coefficients giving the shape and fraction of the YMR population surface brightness profile with respect to the adopted total surface brightness profile from Noyola & Gebhardt (2006); and three Legendre coefficients giving the shape and fraction of the IMR population surface brightness profile. In Figure 10 we show the posterior distributions of the main parameters of this population-dynamical model. It is not possible to show all covariances in one figure, so we show the covariances between parameters pertaining to the Jeans
<figure><img src="image_13.png" /><figcaption>**Figure 10.** MCMC corner plot showing the posterior distributions of the three-population-dynamical model parameters (Section 3.3). Covariances between parameters pertaining to the Jeans equations and isochrone fits are plotted separately.</figcaption></figure>

equations and the isochrone fits separately, although it is a single unified model. The best-fit values of the main model parameters are summarized in Table 2. The posterior of the two Legendre polynomial coefficients, describing the radial fractions of the YMR and IMR populations, are shown in Figure 17 in Appendix A.

We show the best-fit $\Upsilon(r)$ profile for this model in Figure 11. It is very similar to what we found for the two-population model qualitatively and quantitatively (Figure 6). With this model we estimate essentially the same total mass of Sgr's NSC as with the two-population model with radially varying $\Upsilon - 1.60 \pm 0.07 \times 10^6 M_{\odot}$ out to its tidal radius ($\langle\Upsilon\rangle = \frac{M_{\text{tot}}}{L_{\text{tot}}} = 1.96 \pm 0.08 M_{\odot}/L_{\odot}$), of which 65% ($1.04 \pm 0.05 \times 10^6 M_{\odot}$) belongs to the OMP population, 20% ($3.2 \pm 0.2 \times 10^5 M_{\odot}$) belongs to the YMR population, and 15% ($2.4 \pm 0.2 \times 10^5 M_{\odot}$) belongs to the IMR population. Table 5 gives the mass estimates and $\Upsilon_V$ of the three populations locked within several key radial points.

The predictions of this model for the velocity dispersion and rotation profiles of the three different populations are shown in Figure 12. According to the Jeans model, similar to the OMP, the IMR population has a negligible angular momentum, while the YMR population is the only one that owes a significant

fraction of its dynamical support to ordered motions. The latter also has the lowest velocity dispersion. It is interesting to note that there seems to be some evidence of counter rotation in the outer regions of the binned radial velocity profile of the IMR stars, which would make the Sgr dSph another example of a galaxy with a kinematically decoupled core (De Rijcke et al. 2004; Johnston et al. 2018; Fahrion et al. 2019). The discrete Jeans model, however, has not picked up on this, predicting no rotation of the IMR population. This could be due to limitations of the model as the fit is dominated by the denser central parts, where there appears to be no rotation. We caution that due to the overall small number of stars associated with the IMR population, the uncertainties of the associated derived quantities are also considerably higher than for the YMR and OMP populations. In addition, the stars with high probability of belonging to the IMR population have a relatively high velocity dispersion in the outer region of M54 and an overall flatter radial dispersion profile. Our dynamical model does predict a flatter velocity dispersion profile for the IMR stars.

A flatter velocity dispersion profile of the IMR stars could be expected if they also have more extended density distribution and thus, the stars seen at small radii have preferably larger
<table><caption>**Table 5**
Cumulative Dynamical Mass Estimates and $\Upsilon_V$ of M54's Populations According to the Model Presented in Section 3.3</caption><thead><tr><th></th><th>OMP<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>IMR<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>YMR<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>Total<br>(10<sup>6</sup> M<sub>☉</sub>)</th><th>Υ<sub>OMP</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th><th>Υ<sub>IMR</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th><th>Υ<sub>YMR</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th><th>Υ<sub>Tot</sub><br>(M<sub>☉</sub>/L<sub>☉</sub>)</th></tr></thead><tbody><tr><td>r<sub>c</sub></td><td>0.07 ± 0.01</td><td>0.02 ± 0.01</td><td>0.04 ± 0.01</td><td>0.13 ± 0.01</td><td>1.38 ± 0.06</td><td>1.38 ± 0.24</td><td>1.37 ± 0.10</td><td>1.39 ± 0.12</td></tr><tr><td>r<sub>h</sub></td><td>0.56 ± 0.02</td><td>0.09 ± 0.01</td><td>0.22 ± 0.01</td><td>0.87 ± 0.02</td><td>1.65 ± 0.04</td><td>1.51 ± 0.25</td><td>1.51 ± 0.10</td><td>1.59 ± 0.04</td></tr><tr><td>r<sub>fov</sub></td><td>0.90 ± 0.04</td><td>0.15 ± 0.01</td><td>0.27 ± 0.02</td><td>1.32 ± 0.04</td><td>1.92 ± 0.03</td><td>1.78 ± 0.24</td><td>1.65 ± 0.10</td><td>1.83 ± 0.06</td></tr><tr><td>r<sub>t</sub></td><td>1.04 ± 0.05</td><td>0.24 ± 0.02</td><td>0.32 ± 0.02</td><td>1.60 ± 0.07</td><td>2.01 ± 0.08</td><td>2.21 ± 0.17</td><td>1.74 ± 0.10</td><td>1.96 ± 0.08</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** M54 best-fit radial $\Upsilon_V$ profile, according to the three-population-dynamical model described in Section 3.3, resulting from dividing the best-fit mass and luminosity MGE models of the system. The median $\Upsilon_V$ profile is indicated with a thick curve, while the thin curves represent random draws from the posterior distribution. The red dashed curve shows the median $\Upsilon_V$ profile according to the two-population-dynamical model described in Section 3.2. The three black vertical dashed lines correspond to M54's core, half-light, and tidal radii (GC catalog of Harris 1996, 2010 edition), while the red vertical dashed line indicates the radial limit of our kinematic data. The horizontal dashed line corresponds to the average M/L ratio for the entire system, e.g., the total mass divided by total luminosity out to the tidal radius.</figcaption></figure>

physical radii and appear close to the center only in projection. We found this indeed to be the case in Paper I, and here we put this finding to a physical test with the population-dynamical model. As mentioned above, we fit for the radial density profiles of the three-population components within the dynamical model, where their radial luminosity fractions are described by 2nd order Legendre polynomials with the constraint that their sum is equal to the observed surface brightness profile of M54. We show estimated radial density profiles of the three populations in Sgr's NSC and their respective cumulative functions, as well as the population fractions as a function of radius and their cumulative masses in Figure 13. We find that indeed the IMR population has the most extended surface density profile of the three, while the YMR population remains the most centrally concentrated. Our dynamical model predicts that the YMR population is slightly flattened ($q \approx 0.85 \pm 0.04$) as expected from its higher degree of rotation, the OMP one ($q \approx 0.93 \pm 0.04$) appears to be highly spherical with a negligible degree of flattening, and the IMR falls in between ($q \approx 0.89 \pm 0.07$).

In Figure 14 we show how the stars of the three populations are separated probabilistically by our model on the CMD and plot the best-fit isochrones. We estimate the reddening toward M54 to be $E(B-V) = 0.11$ mag. The YMR population has a mean age $2.53 \pm 0.01$ Gyr, the IMR population $5.37 \pm 0.06$ Gyr,

and the OMP one $14.14 \pm 0.13$ Gyr according to the Dartmouth isochrone models. We find that the photometric uncertainties in Sarajedini et al. (2007) catalog and the metallicity spreads appear sufficient to explain the CMD broadening, except for the YMR population, where we need an additional color uncertainty of $\sigma_{V-I}^{YMR} = 0.05$ mag. The median color uncertainty of the photometric catalog is 0.028 mag.

As for the metallicities, we find $\langle[\text{Fe/H}]\rangle_{\text{YMR}} = -0.21 \pm 0.01$ dex, with an intrinsic spread $\sigma_{[\text{Fe/H}]}^{\text{YMR}} = 0.17 \pm 0.01$ dex; $\langle[\text{Fe/H}]\rangle_{\text{IMR}} = -0.47 \pm 0.01$ dex, with an intrinsic spread $\sigma_{[\text{Fe/H}]}^{\text{IMR}} = 0.26 \pm 0.02$ dex; $\langle[\text{Fe/H}]\rangle_{\text{OMP}} = -1.21 \pm 0.01$ dex, with an intrinsic spread $\sigma_{[\text{Fe/H}]}^{\text{OMP}} = 0.24 \pm 0.01$ dex. We find significant intrinsic metallicity spreads in all three populations and while part of these could be physical in nature for the IMR and OMP populations, it is difficult to explain a physical metallicity spread for the YMR population. We note that $\sigma_{[\text{Fe/H}]}^{\text{YMR}}$ is the lowest of the three populations and is likely dominated by systematic factors like underestimated metallicity measurement errors (see Section 5.1 in Husser et al. 2016).

Note also that quoted uncertainties of the best-fit parameters are purely statistical in nature and they are very low due to the large number of stars that are fitted simultaneously. The real uncertainties are dominated by systematic effects in the Dartmouth isochrone models, the photometry, the metallicity estimates, $\alpha$-enhancement, the assumed distance, etc., which exploration is beyond the scope of this work.

Overall, the results of the three-population model in terms of M/L, rotation signatures and dispersion profiles are consistent with our previous model. However, the three-component model can be used to better separate two MR populations in M54, which we find to have different spatial distribution, kinematics, mean ages, and mean metallicities. They likely have different origins, where IMR stars belong to the most central field population of Sgr, while the rapidly rotating YMR population is formed in an in situ burst of star formation in the nucleus of the dwarf galaxy.

# 4. Discussion

## 4.1. Dynamical Model Comparison

We presented three dynamical models in this work, exploring the effects of different gravitational potential parameterization and changing the number of population components. A summary of the estimated total mass and mass contributions of the different populations out to the tidal radius according to the three dynamical models is given in Table 6.

In Sections 3.1 and 3.2, we considered two-population components (MP and MR), but changed the assumptions for the $\Upsilon$ relation—two separate constant M/L ratios for the two populations (Section 3.1) versus a common radially varying M/L ratio (Section 3.2). A problem of the former model is that it attributes almost all of the system's mass to a single population
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Predictions for the rotation (left panel) and the velocity dispersion (right panel) for the three-population-dynamical model (Section 3.3) drawn from the model posterior (YMR—red lines, IMR—green lines, OMP—blue lines), compared to the observed binned profiles (red, green, and blue dots with error bars, respectively), using stars with >50% probability of belonging to either of the three populations. The horizontal error bars indicate the size of each bin. The vertical dashed line in the right panel shows the border between the NFM and WFM MUSE observations.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Top left panel: radial distribution functions of the three-population-dynamical components of M54 (Section 3.3) as drawn from the model posterior (red lines for the YMR, green lines for the IMR, and blue lines for the OMP populations). The solid black line shows the combined surface brightness profile. Top left panel: radial cumulative distribution function of the three-population-dynamical components of M54. Bottom left panel: relative fractions of the three populations as a function of radius. Bottom right panel: cumulative masses of the three populations and the entire system (black lines). The black dashed lines indicate the core, half-light, and tidal radii of the system and the red dashed line indicates the limit of the MUSE field of view.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Top: radial velocity vs. metallicity plot of all M54 member stars in the sample (left panel) and a CMD (right panel), color coded by the probability of the stars belonging to the YMR (red), the IMR (green), and the OMP (blue) populations according to the three-population-dynamical model (Section 3.3). Bottom: separate CMDs of the observed stars, color coded according to the probability to belong to each of the three populations: YMR (shades of red), IMR (shades of green), and OMP (shades of blue). A darker color indicates higher probability. The best-fit Dartmouth isochrones are also plotted together with the HB model.</figcaption></figure>

<table><caption>**Table 6**

Population-dynamical Models Information Criteria, Population Mass Fractions, and Total Mass Estimates at the Tidal Radius</caption><thead><tr><th>Model</th><th>ln P<sub>max</sub></th><th>N<sub>par</sub></th><th>AIC</th><th>BIC</th><th>M<sub>YMR</sub></th><th>M<sub>IMR</sub></th><th>M<sub>OMP</sub></th><th>M<sub>tot</sub><br>10<sup>6</sup> M<sub>☉</sub></th></tr></thead><tbody><tr><td>Section 3.1</td><td>-43110</td><td>24</td><td>86268</td><td>86438</td><td>4%</td><td>...</td><td>96%</td><td>1.26 ± 0.03</td></tr><tr><td>Section 3.2</td><td>-43105</td><td>26</td><td>86262</td><td>86447</td><td>33%</td><td>...</td><td>67%</td><td>1.58 ± 0.07</td></tr><tr><td>Section 3.3</td><td>-29380</td><td>33</td><td>58826</td><td>59060</td><td>20%</td><td>15%</td><td>65%</td><td>1.60 ± 0.07</td></tr></tbody></table>

component. Both models have significantly different predictions for the total mass of M54.

The dynamical model described in Section 3.3 utilizes the gravitational potential parameterization with varying $Υ$ and adds an additional population component. It is fitted to a data set with expanded dimensions because we included photometric information in addition to the metallicity. The resulting total mass estimate for M54 is practically the same as in the two-population-dynamical model with radially varying $Υ$, which did not include photometric data in the fit (Table 6). However, the three-population-dynamical model brings additional insight into the different origins of the MR stars in Sgr's nuclear cluster, which could not be detected with a simpler

model, mainly their different density distributions and angular momentum.

In Figure 15, we compare the best-fit values for the flattening, rotation, and anisotropy of the modeled populations in the three dynamical models. There are no apparent discrepancies between the values of these parameters across the different models. In all three models, MR stars are predicted to rotate faster, have a slightly flatter distribution, and be more tangentially anisotropic than MP stars. The latter, on the other hand, are more radially anisotropic. One could also note that the separation of the YMR and IMR components from the MR population in the three-population model, strengthens the dynamical differences between the most MR and MP stars in
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Comparison between the best-fit values for the flattening (top panel), luminosity weighted rotation in arbitrary units (middle panel), and anisotropy (bottom panel) of M54's modeled populations—YMR (red), IMR (orange), and OMP (blue), according to the three dynamical models described in Section 3.1 (two populations with constant independent $\Upsilon_k$), Section 3.2 (two populations with radially varying $\Upsilon(r)$), and Section 3.3 (three populations with radially varying $\Upsilon(r)$), respectively.</figcaption></figure>

M54. In the three-population model, the YMR component has even higher angular momentum and ellipticity than in the two-population model, where the YMR and IMR components are considered together. Interestingly, in the three-population model, the YMR population is consistent with pure isotropic rotation $\kappa^{\text{YMR}} = -0.98 \pm 0.17$ and $\beta_z^{\text{YMR}} = 0.00 \pm 0.08$ (see Table 2).

In Table 6, we compare the goodness of fit of the three models by looking at the median maximum likelihood ($P_{\text{max}}$) reached by the MCMC fit, the number of free parameters ($N_{\text{free}}$) of each model, and the Akaike and Bayesian information criteria (AIC and BIC). $\text{AIC} = 2 * N_{\text{free}} - 2 * \ln P_{\text{max}}$ is a relative estimator of the quality of the fit, which penalizes the model for the number of free parameters to avoid overfitting. Models with lower AIC values are preferred. Similarly $\text{BIC} = N_{\text{free}} * \ln(n) - 2 * \ln P_{\text{max}}$ is an alternative estimator for the quality of the fit, which also takes into account the sample size ($n$). In our case $n = 8927$ is the number of stars with good quality velocity and metallicity measurements that we use to fit the dynamical models. By its definition, BIC

penalizes a higher number of free parameters stronger than AIC.

It is evident from Table 6 that the three-population-dynamical model formally provides a significantly better goodness of fit to the data, while the two-population-dynamical models with different gravitational potential treatment have comparable performance. Note, however, that the two and three-population model goodness of fit comparison is not entirely fair because the fit was performed on different data sets (including and excluding photometric information).

## 4.2. Comparison with N-body Models

In this section, we compare the results of our Jeans population-dynamical models with the findings by Baumgardt (2017), based on N-body simulations, who find M54's dynamical mass to be $M = 1.62 \pm 0.03 \times 10^6 M_{\odot}$ ($\Upsilon_v = 2.18 \pm 0.20$). This value is in excellent agreement with our derived mass of the entire nuclear system in our three-population model described in Section 3.3 ($M = 1.60 \pm 0.07 \times 10^6 M_{\odot}$; $\Upsilon_v = 1.96 \pm 0.08$) and for the two-population model with a radially varying $M/L$ ratio described in Section 3.2 ($M = 1.58 \pm 0.07 \times 10^6 M_{\odot}$; $\Upsilon_v = 1.92 \pm 0.08$). Some small differences between the two studies, like the lack of systemic rotation in the N-body model, the assumed distance to M54 (23.5 kpc Baumgardt 2017 versus 26.5 kpc; this work) and its luminosity ($8.1 \times 10^5 L_{\odot}$ Baumgardt 2017 versus $8.5 \times 10^5 L_{\odot}$ this work) are enough to explain the small discrepancy in the $M/L$ ratio.

We only find a discrepancy between M54's mass estimates with respect to the N-body analysis for our dynamical model with constant $M/L$ ratios for the two populations (Section 3.1; $M = 1.26 \pm 0.03 \times 10^6 M_{\odot}$; $\Upsilon_v = 1.51 \pm 0.04$), where we find a significantly lower figure. This difference emphasizes the importance for a correct parameterization of the gravitational potential and shows that the radial variation of $\Upsilon$ cannot be ignored.

## 4.3. A GC and a Nucleated Dwarf

In this section, we discuss the results of our Jeans population-dynamical models in the context of historical population studies of Sgr's nuclear cluster. Monaco et al. (2005) present the idea that the Sgr dwarf is a nucleated galaxy with a central overdensity (cusp) of field stars that is coincident, but independent of the presence of a central GC. They use a CMD to separate both stellar structures, which reflect the two dominant components of the Sgr NSC—the MP and MR stars. Additional work by Bellazzini et al. (2008) and Paper II show that the two populations also share a systemic velocity and thus are clearly colocated.

Monaco et al. (2005) show that the surface brightness profile of the entire system can be well described with two King (1962) profiles—the MP GC stars follow a King profile with a core radius $r_c = 6.46$ (0.9 pc) and a tidal radius $r_t = 7.5$ (62 pc) from Trager et al. (1995), while the MR stars are significantly more extended with $r_c = 12.6$ (1.7 pc) and $r_t = 16.7$ (138 pc), and ~30 times less dense in the center than the GC population.

In Figure 16, we compare these two King profiles with our dynamical separation of three stellar structures—the YMR, IMR, and OMP populations. Strictly speaking, our equivalent of the extended population from Monaco et al. (2005) would be the MR population from the two-population-dynamical model
<figure><img src="image_19.png" /><figcaption>**Figure 16.** King profiles of MP stars (blue line Trager et al. 1995) and MR stars (red line Monaco et al. 2005). The combined profile is plotted with a dashed black line and is compared to the surface brightness profile of Noyola & Gebhardt (2006, dashed green line).</figcaption></figure>

(Section 3.2), which is also more spatially extended, but our predicted central density of the MR stars is slightly higher than the central density of the MP stars (Figure 8). However, we have already shown in this work and in Papers I and II, that the MR stars are not a homogeneous population, but rather at least two different stellar structures with different kinematics and origin. In this case, the IMR stars from the three-population model (Section 3.3) are a better representation of the MR nucleus described by Monaco et al. (2005). Figure 16 shows that there is a generally good agreement between the density distributions of the MP stars from Monaco et al. (2005) and our OMP population. The MR Sgr nucleus from Monaco et al. (2005) also follows a very similar density profile to our IMR population, except in the innermost 10", where we find a significantly higher density. The discrepancy could be due to lack of sufficient spatial resolution in the work by Monaco et al. (2005). We find the central density of IMR stars to be five to six times lower than the central density of OMP stars.

## 4.4. The Extended SFH of Sgr

The Sgr dSph has a very extended SFH (de Boer et al. 2015) and thus we expect that the field population also includes stars that have essentially the same ages and metallicities as the OMP NSC population. However, in the population-dynamical models presented here, we model the central overdensity of Sgr field stars (the IMR component) as a single stellar population with a mean age and metallicity. We do fit for its intrinsic metallicity spread and we find $\sigma_{[Fe/H]}^{\text{IMR}} = 0.26 \pm 0.02$ dex, the largest of the three populations. The intrinsic metallicity spread of the OMP population is also significantly larger than expected for a GC ($\sigma_{[Fe/H]}^{\text{OMP}} = 0.25 \pm 0.01$ dex). There are different possible explanations for the increased metallicity spread of the OMP stars—either this sample is contaminated by MP Sgr field stars, or it is the post-merger remnant of multiple GCs with slightly different metallicities, or both. In addition, part of the detected metallicity spread among the individual stellar populations could also be due to underestimated [Fe/H] measurement uncertainties, as discussed in Husser et al. (2016).

In Paper I, we measured individual stellar ages and showed that IMR stars have the largest intrinsic age spread ($1.16 \pm 0.07$ Gyr), although this figure still does not include

any MP stars. Here, we chose to work with the more basic CMD information (magnitudes and colors), instead of relying on individual age estimates, which prevents us from deriving directly age spread estimates for the model-defined populations in M54. However, we do not detect any additional broadening of the CMD than expected from the photometric uncertainties and metallicity spreads for either the OMP or IMR population. We do not expect to be able to distinguish between photometric broadening on the CMD caused by intrinsic age or metallicity spreads, as the effects are highly degenerate, especially at older ages.

Both the IMR and OMP populations appear to have very similar kinematic properties in our dynamical model (i.e., isotropic and nonrotating); hence, we expect that we cannot discriminate between MP field and GC stars kinematically. If these populations also overlap in age and metallicity space, the only discriminating factor remains their different radial density distributions. The Sgr field is indeed more radially extended than the GC population, considering the radial distribution of its MR tail (the IMR population); however, the OMP stars are significantly more abundant than the IMR stars at all radii in our data coverage, which makes it impossible to separate them in a meaningful way.

We thank Holger Baumgardt for valuable discussions about the shape of M54's gravitational potential. We thank Tim-Oliver Husser for providing the MUSE LSF corrected PHOENIX spectral library. This research has made use of NASA's Astrophysics Data System. G.v.d.V. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 724857 (Consolidator Grant ArcheoDyn). A.M.B. acknowledges funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 895174.

# 5. Summary

This is the third in a series of papers aiming to reexamine the massive star cluster M54 as Sgr's dwarf galaxy's nuclear cluster. In the first two papers of the series we published results based on the extensive MUSE seeing-limited WFM mosaic data, obtained by our team during observing run 095.B-0585 (A) (P.I.: Lützgendorf). Our analysis then included $> 7400$ individual stellar spectra with a signal-to-noise ratio (S/N) $> 10$ $px^{-1}$ extracted from the MUSE cubes with pampelmuse (Kamann et al. 2013). In this third paper of the series we supplement the original stellar sample with additional WFM-AO and NFM-AO observations in the center of the system, taken during the science verification of the AO capabilities of MUSE (P.I.: Alfaro-Cuello) for a total of $\sim 9000$ unique stellar targets with S/N $> 10$ $px^{-1}$ and improved spatial resolution in the inner region of M54.

In Paper I, we obtained spectroscopic metallicities and measured individual stellar ages from isochrone fits for the majority of the stars with extracted spectra. We classified them into three distinct stellar populations that we called OMP, IMR, and YMR and showed that they have different spatial distributions in the core of Sgr. In Paper II, we explored the kinematic properties of these three stellar populations and showed that they exhibit different rotation and velocity dispersion profiles. We offered a formation scenario for Sgr's
NSC, which consists of stars with different origins and hence different population and kinematic characteristics:

i) A mixture of globular cluster stars (OMP) that belong to at least one massive, but possibly a merger of several GCs, due to their large metallicity spread

ii) Stars formed recently in situ (YMR) that are very centrally concentrated and still have a high degree of rotation;

iii) Stars that belong to the inner Sgr field population (IMR) that have a large age spread, very extended spatial distribution, and a relatively flat velocity dispersion profile.

In this work we explore the dynamical imprints of this formation scenario and the interplay between the different populations in a common gravitational potential. We analyzed all individual MUSE stellar spectra with SPEXXY, utilizing the phoenix stellar library to derive radial velocities and metallicities. This includes a reanalysis of the old sample with the slightly different methodology adopted in the current paper to get a truly uniform sample. We use a dynamical modeling approach, based on the Jeans equations. Each of the dynamical components in the model is characterized by its population properties (mean metallicity and age), spatial distribution, and velocity moments (rotation and velocity dispersion profiles), which we fit for simultaneously. Our population-dynamical models are successful in estimating the joint probability of each star in the spectroscopic sample of belonging to either one of the dynamical components of the model, or to the Milky Way foreground, based on its observed quantities (radial velocity metallicity, photometry, and their respective uncertainties).

Overall, we show that simple population-dynamical models, based on the Jeans equations, can explain the majority of observed properties of M54's complex stellar populations simultaneously and self-consistently. We confirm our previous findings and especially emphasize that the MR stars in M54 have a heterogeneous origin. The YMR population must have formed from in situ star formation in a gaseous disk. It is the most centrally concentrated, slightly more flattened than the rest, and with the highest angular momentum. The IMR population corresponds to Sgr dSph's inner field population. It is very extended spatially, follows a relatively flat velocity dispersion profile, and has a considerable metallicity spread. Finally, the OMP population has all the characteristics that are typical for globular clusters. Such a mixed nuclear formation mechanism has been described in the literature (Neumayer et al. 2011; den Brok et al. 2014b; Antonini et al. 2015; Cole et al. 2017; Neumayer et al. 2020; Fahrion et al. 2022).

We also look into M54's mass and *M/L* ratio (Υ) radial profiles. We conclude that taking into account the radial variations of the *M/L* ratio is important to correctly reproduce the gravitational potential realistically. We find that Υ<sub>V</sub> is U-shaped with a minimum around the core radius, followed by a monotonic increase outward (Figure 11). Multiple factors contribute to the complex Υ profile. Mass-segregated dark remnants and the centrally concentrated young stars have opposing effects on the *M/L* ratio in the inner regions of the NSC, while the strong increase outward could be due to an increased number of low mass stars and non-negligible dark matter contribution. The total dynamical mass enclosed within the tidal radius (76 pc) of the entire system and its different stellar components are summarized in Table 6. Our results are in excellent agreement with the *N*-body simulation studies by Baumgardt (2017).

In a forthcoming work (M. Alfaro-Cuello et al., in preparation) we will address the question whether an intermediate-age black hole exists at the core of M54, based on the dynamical models presented here and the work by Aros et al. (2020). These results provide a benchmark for studying the disruption of nucleated satellite galaxies and the formation of ultracompact dwarfs.

We thank Holger Baumgardt for valuable discussions about the shape of M54's gravitational potential. We thank Tim-Oliver Husser for providing the MUSE LSF corrected PHOENIX spectral library. This research has made use of NASA's Astrophysics Data System. G.v.d.V. acknowledges 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 Archeo Dyn). A.M.B. acknowledges funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 895174.

**Appendix A**

**Corner Plots of the Individual Population Surface Brightness Distribution Parameterizations**

Figure 17 shows the covariances of the 10 additional para-meters describing the MGE fractions contributing to the MP population in the two-population model (Section 3.1) and the Legendre coefficients parameterizing the MGE fractions in the three-population model (Section 3.3).
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Top: corner plot of the individual fractions of the MGE components contributing to the MP population according to the best-fit two-population-dynamical model presented in Section 3.1. Bottom: corner plot of the Legendre polynomial coefficients used to compute the YMR (**left**) and IMR (**right**) fractions of the surface brightness MGE in the three-population-dynamical model presented in Section 3.3.</figcaption></figure>
## Appendix B
Complete M54 Data Set

The complete M54 data set is available in a machine-readable format. Table 7 provides a short sample of the complete version to guide its reuse.

<table><caption>**Table 7**
Complete M54 Data Set Used in This Work</caption><thead><tr><th>ID</th><th>R.A.<br>(deg)</th><th>Decl.<br>(deg)</th><th>V<sub>r</sub><br>(km s<sup>-1</sup>)</th><th>Err. V<sub>r</sub><br>(km s<sup>-1</sup>)</th><th>[Fe/H]</th><th>Err. [Fe/H]</th><th>F606W<br>(mag)</th><th>F814W<br>(mag)</th><th>Err. F606W<br>(mag)</th><th>Err. F814W<br>(mag)</th><th>Pointing</th></tr></thead><tbody><tr><td>36954</td><td>283.7653198</td><td>-30.4947300</td><td>141.419</td><td>0.600</td><td>-1.043</td><td>0.017</td><td>17.735</td><td>16.882</td><td>0.005</td><td>0.012</td><td>WFM-noAO</td></tr><tr><td>41660</td><td>283.7614136</td><td>-30.4950466</td><td>130.990</td><td>2.332</td><td>-1.251</td><td>0.059</td><td>19.202</td><td>18.414</td><td>0.006</td><td>0.016</td><td>WFM-noAO</td></tr><tr><td>40701</td><td>283.7621765</td><td>-30.4926395</td><td>153.876</td><td>1.958</td><td>-1.360</td><td>0.097</td><td>18.286</td><td>18.075</td><td>0.032</td><td>0.015</td><td>WFM-noAO</td></tr></tbody></table>

**Note.** Here we show only the first three entries, the full table is available online. The stellar IDs, coordinates, and photometry are from the HST catalog of Sarajedini et al. (2007). The radial velocities and metallicities are from this work.
(This table is available in its entirety in machine-readable form.)

## ORCID iDs

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

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

Nora Lützgendorf https://orcid.org/0000-0002-4034-0080

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

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

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

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

Karina T. Voggel https://orcid.org/0000-0001-6215-0950

Iskren Y. Georgiev https://orcid.org/0000-0001-8471-6679

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

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

## References

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

Alfaro-Cuello, M., Kacharov, N., Neumayer, N., et al. 2020, ApJ, 892, 20

Antonini, F., Barausse, E., & Silk, J. 2015, ApJ, 812, 72

Antonini, F., Capuzzo-Dolcetta, R., Mastrobuono-Battisti, A., & Merritt, D. 2012, ApJ, 750, 111

Aros, F. I., Sippel, A. C., Mastrobuono-Battisti, A., et al. 2020, MNRAS, 499, 4646

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

Bassino, L. P., Muzzio, J. C., & Rabolli, M. 1994, ApJ, 431, 634

Baumgardt, H. 2017, MNRAS, 464, 2174

Bekki, K. 2015, ApJL, 812, L14

Bekki, K., Couch, W. J., Drinkwater, M. J., & Shioya, Y. 2003, MNRAS, 344, 399

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

Bianchini, P., Sills, A., & Miholics, M. 2017a, MNRAS, 471, 1181

Bianchini, P., Sills, A., van de Ven, G., & Sippel, A. C. 2017b, MNRAS, 469, 4359

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

Cappellari, M. 2002, MNRAS, 333, 400

Cappellari, M. 2008, MNRAS, 390, 71

Cappellari, M. 2012, arXiv:1211.7009

Capuzzo-Dolcetta, R., & Miocchi, P. 2008a, MNRAS, 388, L69

Capuzzo-Dolcetta, R., & Miocchi, P. 2008b, ApJ, 681, 1136

Carlberg, R. G., & Grillmair, C. J. 2022, ApJ, 935, 14

Carretta, E., Bragaglia, A., Gratton, R., D'Orazi, V., & Lucatello, S. 2009, A&A, 508, 695

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

Cole, D. R., Debattista, V. P., Varri, A.-L., Adam, M., & Seth, A. C. 2017, MNRAS, 466, 2895

de Boer, T. J. L., Belokurov, V., & Koposov, S. 2015, MNRAS, 451, 3489

De Rijcke, S., Dejonghe, H., Zeilinger, W. W., & Hau, G. K. T. 2004, A&A, 426, 53

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

den Brok, M., van de Ven, G., van den Bosch, R., & Watkins, L. 2014a, MNRAS, 438, 487

Di Cintio, A., Mostoghiu, R., Knebe, A., & Navarro, J. 2021, MNRAS, 506, 531

Dotter, A., Chaboyer, B., Jevremović, D., et al. 2008, ApJS, 178, 89

Emsellem, E., Monnet, G., & Bacon, R. 1994, A&A, 285, 723

Fahrion, K., Leaman, R., Lyubenova, M., & van de Ven, G. 2022, A&A, 658, A172

Fahrion, K., Lyubenova, M., van de Ven, G., et al. 2019, A&A, 628, A92

Fahrion, K., Lyubenova, M., van de Ven, G., et al. 2021, A&A, 650, A137

Fétick, R. J. L., Fusco, T., Neichel, B., et al. 2019, A&A, 628, A99

Filippenko, A. V., & Ho, L. C. 2003, ApJL, 588, L13

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

Gnedin, O. Y., Ostriker, J. P., & Tremaine, S. 2014, ApJ, 785, 71

Goodman, J., & Weare, J. 2010, CAMCS, 5, 65

Graham, A. W., & Spitler, L. R. 2009, MNRAS, 397, 2148

Hannah, C. H., Seth, A. C., Nguyen, D. D., et al. 2021, AJ, 162, 281

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

Hopkins, P. F., Murray, N., Quataert, E., & Thompson, T. A. 2010, MNRAS, 401, L19

Husser, T.-O., Kamann, S., Dreizler, S., et al. 2016, A&A, 588, A148

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

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

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

Jeans, J. H. 1915, MNRAS, 76, 70

Johnston, E. J., Hau, G. K. T., Coccato, L., & Herrera, C. 2018, MNRAS, 480, 3215

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

Kamann, S., Dalessandro, E., Bastian, N., et al. 2020, MNRAS, 492, 966

Kamann, S., Wisotzki, L., & Roth, M. M. 2013, A&A, 549, A71

King, I. 1962, AJ, 67, 471

Layden, A. C., & Sarajedini, A. 2000, AJ, 119, 1760

Leibundgut, B., Bacon, R., Bian, F., et al. 2019, Msngr, 176, 16

Lotz, J. M., Telford, R., Ferguson, H. C., et al. 2001, ApJ, 552, 572

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

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

McLaughlin, D. E., & van der Marel, R. P. 2005, ApJS, 161, 304

Mihos, J. C., & Hernquist, L. 1994, ApJL, 437, L47

Milosavljević, M. 2004, ApJL, 605, L13

Monaco, L., Bellazzini, M., Ferraro, F. R., & Pancino, E. 2005, MNRAS, 356, 1396

Neumayer, N., Seth, A., & Böker, T. 2020, A&ARv, 28, 4

Neumayer, N., & Walcher, C. J. 2012, AdAst, 2012, 709038

Neumayer, N., Walcher, C. J., Andersen, D., et al. 2011, MNRAS, 413, 1875

Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2019, ApJ, 872, 104

Norris, M. A., Escudero, C. G., Faifer, F. R., et al. 2015, MNRAS, 451, 3615

Noyola, E., & Gebhardt, K. 2006, AJ, 132, 447

Oh, K. S., & Lin, D. N. C. 2000, ApJ, 543, 620
Percival, S. M., Salaris, M., Cassisi, S., & Pietrinferni, A. 2009, ApJ, 690, 427

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

Pfeffer, J., Lardo, C., Bastian, N., Saracino, S., & Kamann, S. 2021, MNRAS, 500, 2514

Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523

Rossa, J., van der Marel, R. P., Böker, T., et al. 2006, AJ, 132, 1074

Sarajedini, A., Bedin, L. R., Chaboyer, B., et al. 2007, AJ, 133, 1658

Schinnerer, E., Böker, T., Meier, D. S., & Calzetti, D. 2008, ApJL, 684, L21

Seth, A., Agüeros, M., Lee, D., & Basu-Zych, A. 2008, ApJ, 678, 116

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

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

Trager, S. C., King, I. R., & Djorgovski, S. 1995, AJ, 109, 218

Tremaine, S. D., Ostriker, J. P., & Spitzer, L., Jr. 1975, ApJ, 196, 407

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

Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978

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

Walcher, C. J., Böker, T., Charlot, S., et al. 2006, ApJ, 649, 692

Watkins, L. L., van de Ven, G., den Brok, M., & van den Bosch, R. C. E. 2013, MNRAS, 436, 2598

Wrobel, J. M., Greene, J. E., & Ho, L. C. 2011, AJ, 142, 113

Zhu, L., van de Ven, G., Watkins, L. L., & Posti, L. 2016, MNRAS, 463, 1117