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

# A Deep View into the Nucleus of the Sagittarius Dwarf Spheroidal Galaxy with MUSE.
## I. Data and Stellar Population Characterization

M. Alfaro-Cuello<fnref n="1" />, N. Kacharov<fnref n="1" />, N. Neumayer<fnref n="1" />, N. Lützgendorf<fnref n="2" />, A. C. Seth<fnref n="3" />, T. Böker<fnref n="2" />, S. Kamann<fnref n="4" />, R. Leaman<fnref n="1" />,
G. van de Ven<fnref n="5" />, P. Bianchini<fnref n="6" />, L. L. Watkins<fnref n="7" />,<fnref n="8" />, and M. Lyubenova<fnref n="8" />

Received 2019 January 23; revised 2019 April 16; accepted 2019 April 18; published 2019 November 19

## Abstract

The center of the Sagittarius dwarf spheroidal galaxy (Sgr dSph) hosts a nuclear star cluster (NSC), M54, which is the only galaxy nucleus that can be resolved into individual stars at optical wavelengths. It is thus a key target for understanding the formation of NSCs and their relation to globular clusters (GCs). We present a large Multi-Unit Spectroscopic Explorer data set that covers M54 out to ~2.5 half-light radius, from which we extracted the spectra of ~6600 cluster member stars. We use these data in combination with *Hubble Space Telescope* photometry to derive age and metallicity for each star. The stellar populations show a well-defined age–metallicity relation, implying an extended formation history for the central region of Sgr dSph. We classify these populations into three groups, all with the same systemic velocity: young metal-rich (YMR; 2.2 Gyr, [Fe/H] = −0.04); intermediate-age metal-rich (IMR; 4.3 Gyr, [Fe/H] = −0.29); and old metal-poor (OMP; 12.2 Gyr, [Fe/H] = −1.41). The YMR and OMP populations are more centrally concentrated than the IMR population, which are likely stars of the Sgr dSph. We suggest that the OMP population is the result of accretion and merging of two or more old and metal-poor GCs dragged to the center by dynamical friction. The YMR is consistent with being formed by in situ star formation in the nucleus. The ages of the YMR population suggest that it may have been triggered into forming when the Sgr dSph began losing its gas during the most recent interaction with the Milky Way, ~3 Gyr ago.

**Key words:** galaxies: dwarf – galaxies: individual (Sgr dSph) – galaxies: nuclei – galaxies: star clusters: individual (M54) – globular clusters: individual (M54)

## 1. Introduction

Nuclear star clusters (NSCs) are the densest stellar systems in the universe (Walcher et al. 2005; Misgeld & Hilker 2011; Norris et al. 2014), with masses from $10^6$ to $10^8 M_\odot$ and half-light radii of about 1–10 pc (Georgiev & Böker 2014; Georgiev et al. 2016). High-resolution *Hubble Space Telescope* (HST) imaging of nearby galaxies has revealed that NSCs are found in $\ge 70\%$ of galaxies across the Hubble sequence (Phillips et al. 1996; Carollo et al. 1998; Böker et al. 2002, 2004; Côté et al. 2006; Turner et al. 2012; Georgiev & Böker 2014), covering a wide range of galaxy stellar masses with an occupation fraction peaking at galaxy masses in the range of $10^8-10^{10} M_\odot$ (Georgiev et al. 2009; Ordenes-Briceño et al. 2018; Sánchez-Janssen et al. 2019). In addition, NSCs are known to coexist with supermassive black holes (SMBHs; Seth et al. 2008; Neumayer & Walcher 2012), as observed in the Milky Way (Genzel et al. 2010; Feldmeier et al. 2014; Schödel et al. 2014a, 2014b).

NSCs may also be relevant for the origin of the most massive globular clusters (GCs), which is still a debated question. Many of the most massive clusters in the Milky Way and M31 show significant spreads in metallicity (Meylan et al. 2001; Fuentes-Carrera et al. 2008; Willman & Strader 2012; Bailin 2018), and the most massive Milky Way GC, $\omega$ Cen ($2.5 \times 10^6 M_\odot$; van de Ven et al. 2006), seems to have age spreads as well (Hilker et al. 2004; Villanova et al. 2014). Terzan 5 is another multicomponent

cluster that resides in the bulge of the Milky Way. It displays a high spread in iron (Massari et al. 2014) and age (Ferraro et al. 2016; Origlia et al. 2019), with a chemical composition similar to the Galactic bulge. This led to the suggestion that they are former NSCs of dwarf galaxies that have been accreted (Böker 2008; Da Costa 2016), as NSCs are known to have extended star formation histories (SFHs; e.g., Walcher et al. 2005; Kacharov et al. 2018). Metallicity spreads alone could also be due to self-enrichment during formation (Bailin 2018); however, age spreads are likely a strong indicator of a stripped NSC. Pfeffer et al. (2014) predict, via cosmological simulations of the Milky Way, that between one and three GCs with masses higher than $10^5 M_\odot$ are tidally striped nuclei, with a high likelihood of remnants above $10^6 M_\odot$.

The most massive compact stellar systems are often referred to as ultracompact dwarfs (UCDs), which blend smoothly in mass and size with the massive GCs, with no clear dividing line between them. For massive UCDs, above $10^7 M_\odot$, a large fraction do appear to be stripped nuclei, as they host SMBHs that make up a large fraction of their mass (Seth et al. 2014; Ahn et al. 2017, 2018; Afanasiev et al. 2018) or have an extended SFH (Norris et al. 2015). At masses below $10^7 M_\odot$, no equivalent evidence has been found outside of $\omega$ Cen’s extended SFH, but stripped NSCs at these masses are expected to be quite common (Georgiev & Böker 2014; Pfeffer et al. 2014).

In this context, the following questions arise: How many of the massive GCs in our Milky Way are former nuclei of low-

<footnote n="1">¹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="2">²European Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="3">³Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA</footnote>

<footnote n="4">⁴Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK</footnote>

<footnote n="5">⁵Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17, A-1180 Wien, Austria</footnote>

<footnote n="6">⁶Observatoire astronomique de Strasbourg, CNRS, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="7">⁷Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="8">⁸European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany</footnote>
mass galaxies? What are the detailed properties (morphology, stellar ages, metallicity, and kinematics) we would expect these systems to have?

To address these questions, we focus here on an ideal target in the nucleated Sagittarius dwarf spheroidal galaxy (Ibata et al. 1994)—hereafter Sgr dSph. This galaxy is currently being disrupted by the tidal field of the Milky Way (Ibata et al. 1997), leaving visible streams as evidence. The original luminosity of the galaxy has been estimated to be ~$10^8 L_\odot$ ($M_v \sim -15$; Niederste-Ostholt et al. 2012).

The NSC of this galaxy was discovered as the GC M54 (NGC 6715) long before the remainder of the galaxy was discovered. M54 has a mass of $(1.41 \pm 0.02) \times 10^6 M_\odot$ (Baumgardt & Hilker 2018), making it the second most massive GC in the galaxy after $\omega$ Cen. It is located in the densest region of the Sagittarius stream at the photometric center of the galaxy (Ibata et al. 1994; Mucciarelli et al. 2017) at a distance of 28.4 kpc (Siegel et al. 2011).

The Sgr dSph NSC is composed of at least two distinct populations, a metal-poor one and a metal-rich one, with the metal-rich population having stars as young as ~2 Gyr, while the metal-poor population is consistent with an age of ~13 Gyr (Monaco et al. 2005a; Siegel et al. 2007; Bellazzini et al. 2008). Some authors have referred to these two components as separate objects, a massive metal-poor GC (M54) and a metal-rich component associated with the galaxy (Sgr, or Sgr NS). However, in the context of extragalactic NSCs (where these two populations could not easily be separated), and given the populations' identical photocenters (Monaco et al. 2005b) and line-of-sight (LOS) velocities (Bellazzini et al. 2008), we prefer to think of these components as two (or more) subpopulations of the Sgr dSph NSC. Hence, we will refer throughout the paper to the different subpopulations in the NSC of the Sgr dSph, considering their stars physical parameters (e.g., metallicity, age).

Monaco et al. (2005b) found that the metal-rich subpopulation presents a cusp with a peak located in an indistinguishable position from the center of the metal-poor population, thus displaying a different radial profile in the inner parts in comparison with the rest of the Sgr dSph (Majewski et al. 2003). Ibata et al. (2009) reported that the metal-poor subpopulation shows a stellar density cusp, together with a peak in the velocity dispersion. The density profiles from both subpopulations are found to be different (Monaco et al. 2005b; Bellazzini et al. 2008).

Siegel et al. (2007), using HST photometry, suggested that at least four discrete stellar populations are present in the nucleus of the Sgr dSph, revealing a very complex SFH. These populations are as follows:

(i) [Fe/H] = -1.8, [α/Fe] = +0.2 with ages of 13 Gyr;

(ii) [Fe/H] = -0.6, [α/Fe] = -0.2 with ages of 4-6 Gyr;

(iii) [Fe/H] = -0.1, [α/Fe] = -0.2 with ages of 2.3 Gyr; and

(iv) [Fe/H] = +0.6, [α/Fe] = 0.0 with ages of 0.1-0.8 Gyr.

Previous spectroscopic studies have focused on the division of the cluster into metal-rich and metal-poor subpopulations (i.e., population (i) above versus all others). Bellazzini et al. (2008) used ~400 stars to show that the systemic velocities of the stars in both metallicity regimes coincide within $\simeq \pm 1.0$ km s$^{-1}$, consistent with previous findings based on samples with a considerably lower number of stars (Da Costa &

Armandroff 1995; Ibata et al. 1997; Monaco et al. 2005a). Despite this coincidence in radial velocities, the velocity dispersion profiles differ significantly for these populations, with the metal-rich stars having a much flatter dispersion profile than the metal-poor stars (Bellazzini et al. 2008). To explain the spatial coincidence and the differing dispersion and surface brightness profiles of the different populations, Monaco et al. (2005b) and Bellazzini et al. (2008) suggested the following possible scenarios: (1) Both were formed in situ at the bottom of the potential well of the galaxy, first the metal-poor stars, followed by the metal-rich stars in subsequent star-forming episodes from enriched gas. (2) The metal-poor stars were an ordinary GC that was driven by dynamical friction to the center of the galaxy, where a nucleus either had already formed independently or formed subsequently.

Using N-body simulations, Bellazzini et al. (2008) showed that the latter scenario is feasible even with a preexisting NSC. More specifically, the dynamical friction inspiral timescale of the metal-poor progenitor GC is <3 Gyr for a wide range of starting radii and initial relative velocities (unlike the other GCs in the Sgr dSph). Furthermore, the end state of these simulations results in systemic velocity differences between the preexisting NSC and migrated GC of <2 km s$^{-1}$.

Studies on the abundances of the metal-poor and metal-rich population and radial variations of the metallicity have also been published. The metal-poor population displays a large spread in the iron content of its stars ($\sigma_{[Fe/H]}$ = 0.186; Carretta et al. 2010b; Willman & Strader 2012). In addition, Carretta et al. (2010b) reported that a sodium-oxygen anticorrelation similar to those seen in other GCs is present in the metal-poor population. On the other hand, the authors do not observe any signatures of this anticorrelation in the metal-rich population. Mucciarelli et al. (2017) found a metallicity gradient in the metal-rich population, with this component becoming more metal-poor at larger radii. They also estimated that the chemical distribution they found for the Sgr dSph is consistent with a progenitor of $M_{DM} = 6 \times 10^{10} M_\odot$, in agreement with the estimate by Gibbons et al. (2017).

Finally, kinematic observations of the metal-poor population have also been used to determine whether it hosts a central massive black hole. Two studies have reported a possible detection based on dynamical modeling ($M_{BH} \sim 10^4 M_\odot$; Ibata et al. 2009; Baumgardt 2017).

All of these characteristics make the nucleus of the Sgr dSph a key object for understanding NSCs. Its proximity enables us to estimate parameters for individual stars, providing valuable information on the nucleus itself and its coevolution with the progenitor galaxy. This object offers the opportunity to better understand the link between the most massive GCs and NSCs as stripped nuclei from low-mass galaxies. The Multi-unit Spectroscopic Explorer (MUSE) spectrograph is an excellent instrument for studying individual stars in GCs (e.g., Husser et al. 2016; Kamann et al. 2018). In this work, we present a powerful MUSE data set centered on M54. With these data we perform a spectroscopic study of over 6500 stars, an order-of-magnitude increase over previous spectroscopic samples.

In this paper we present our data in Section 2, including the extraction of the single stellar spectra. Section 3 presents our analysis, including radial LOS, membership probability, metallicity, and age measurements. Our results are presented in
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Color image obtained from the MUSE data of the Sgr dSph NSC, using synthetic *i*, *r*, and *z* filters. The image is a 4 × 4 mosaic of 16 MUSE pointings, covering ~2.5$R_{\text{HL}}$ ($R_{\text{HL}} = 0.82$; Harris 1996, 2010 edition) of this NSC. The total covered FOV—considering overlaps—is ~3.5 × 3.5, corresponding to 29 pc × 29 pc (1' = 8.27 pc). North is up and east is to the left.</figcaption></figure>

Section 4; we discuss these in Section 5 and conclude in Section 6.

## 2. Data

### 2.1. *Observations and Data Reduction*

The data set was acquired with MUSE (Bacon et al. 2014), an integral field spectrograph mounted at the UT4 of the Very Large Telescope at the Paranal Observatory in Chile, between 2015 June 29 and September 19 in run 095.B-0585(A) (PI: Lützgendorf).

A total of 16 pointings, with a field of view (FOV) of 59.9 × 60.0 each, create a 4 × 4 mosaic that covers a contiguous area out to ~2.5 times the half-light radius of M54 ($R_{\text{HL}} = 0.82$, Harris 1996, 2010 edition), equivalent to ~25 times its core radius ($R_c = 0.09$, Harris 1996, 2010 edition). The estimated tidal radius of M54 is $R_t = 9.87$ (assuming a King-model central concentration of $c = 2.04$; Harris 1996, 2010 edition).

Due to overlaps between neighboring pointings, the observations cover an FOV of ~3.5 × 3.5 (1' corresponds to 8.27 pc at an assumed distance of 28.4 kpc; Siegel et al. 2011). This results

in a total FOV of 29 pc × 29 pc. The data have a wavelength coverage of 4800–9300 Å and a spectral resolution increasing with wavelength from $R \sim 1750$ to $R \sim 3750$. The spatial sampling is 0.2 pixel$^{-1}$. Each field was observed with three exposures, applying offsets in rotation of 90° between them (no dithering) to increase the homogeneity of the data across the FOV. The exposure time for each of the three exposures was chosen to be 400 s for the inner 4 pointings and 610 s for the outer 12 pointings. This was done to avoid saturation of the innermost crowded region with high surface brightness. During the observations, the air mass varied from 1.0 to 1.9, and the seeing between 0.5 and 2. Figure 1 shows the color image of the 4 × 4 mosaic obtained from the MUSE data using synthetic *i*, *r*, and *z* filters.

The data reduction was performed using the MUSE pipeline (version 1.2.1; Weilbacher et al. 2014). The pipeline includes all tasks for the data reduction process: bias creation/ subtraction, flat-fielding, illumination correction, and wavelength calibration. Before combining the single exposures, they were flux-calibrated using a standard star and corrected for the barycentric motion of Earth.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Three examples of the fitting performed by ULySS (Koleva et al. 2009) for different types of stars with S/Ns of 14 (close to minimum of 10), 48, and 100, in the wavelength range of 4701–6800 Å. Top panels: the best fit is represented in blue, data not considered for the fit in red (both panels), and the multiplicative polynomial in turquoise. Bottom panels: residual of the spectrum fitting (black). Overplotted in dashed and solid green are the mean and the $1\sigma$ deviation, respectively.</figcaption></figure>

## 2.2. Stellar Spectra Extraction

We extract individual stellar spectra with PampelMuse (Kamann et al. 2013). This tool takes a photometric catalog as a reference for the positions and magnitudes of the stars in the FOV. It models the change in the point-spread function

(PSF) as a function of wavelength and allows deblending of sources efficiently even in crowded and dense regions. A reference stellar catalog for M54 was published by Siegel et al. (2007) as part of the ACS Survey of Galactic Globular Clusters (Sarajedini et al. 2007). The catalog includes source positions
and F606W (V) and F814W (I) magnitudes with their
respective uncertainties. Its total coverage is larger than the
MUSE data set FOV.

We extract the spectra of ~55,000 stars from the entire field.
For our subsequent analysis we only consider the spectra with
signal-to-noise ratio (S/N) ≥ 10, which are labeled by
PampelMuse as good extracted spectra. Within the entire
MUSE FOV, we extract good spectra for a total of ~7000
different stars in the magnitude range of 22 < F814W < 13
(I band).

In our stellar sample we do not include main-sequence (MS)
stars below the oldest turnoff point since they are fainter than
I ~ 21 and their spectra are typically below our S/N = 10
threshold.

**3. Analysis**

To measure physical stellar parameters, we use University of
Lyon Spectroscopic Analysis Software (ULySS; Koleva et al.
2009), a tool for determination of atmospheric parameters of
stars through interpolating and fitting stellar library templates to
the observed spectrum in the wavelength range of
3900–6800 Å. We note that the spectral wavelength range of
MUSE allows us to fit from 4800 to 6800 Å only, and we
masked the Na I D line region between 5850 and 5950 Å of the
observed spectrum. The observed spectrum is compared with a
library of synthetic spectra with varying metallicity, surface
gravity, and temperature. The best fit is achieved by
interpolation within the synthetic library and χ² minimization.
The fit performed by UlySS is a Levenberg–Marquardt local
minimization for nonlinear parameters, while for the linear
ones it uses a bounded-values least square. The synthetic
spectroscopic grid available in ULySS is built on the basis of
the ELODIE 3.2 library (Wu et al. 2011). It is limited to a
metallicity range of −2.5 < [Fe/H] < 0.5. ULySS also
simultaneously estimates the LOS velocity shifts of the input
spectra. Finally, ULySS derives the uncertainties of the
different physical parameters via Monte Carlo simulations.
For stars in overlapping regions, where multiple spectra are
extracted, the final values for the LOS velocity and metallicity
are the weighted mean and its uncertainty of the multiple
ULySS measurements.

Figure 2 shows three examples of the fitting performed by ULySS for stars with a spectrum of S/Ns of 14 (close to the minimum value of 10), 48, and 100. We use a multiplicative polynomial of order 6, which provides sufficiently good normalization of the observed spectra without the danger of overfitting.

3.1. *Line-of-sight Velocity*

From the spectral fitting performed by ULySS we obtain the
LOS velocities for each star.

Since the space motion of the nucleus can mimic the effect of
an additional rotation component, we correct the LOS velocity
values for this effect, known as perspective rotation, using
Equation (6) in van de Ven et al. (2006). We use the absolute
proper motion of the cluster M54 $V_W = -2.82 \pm 0.11$ mas yr$^{-1}$
and $V_N = -1.51 \pm 0.14$ mas yr$^{-1}$ (Sohn et al. 2015). These
values are in good agreement with the recent values based on
Gaia data by Vasiliev (2019). We obtain correction values
between -0.3 and 0.3 km s$^{-1}$.

<figure><img src="image_4.png" /><figcaption>Figure 3. Top: distribution of the LOS velocities for the extracted stars with S/N ≥ 10 from the entire MUSE FOV. Bottom: radius vs. velocity for our stellar sample color-coded by the cluster membership probability ($P_{\text{Membership}}$). We consider as member stars of the Sgr dSph NSC those with probability ≥70%.</figcaption></figure>

In Appendix A, we discuss the agreement between velocity
measurements derived for stars observed in multiple MUSE
pointings.

3.2. Sgr dSph NSC Membership

The position of the Sgr dSph with respect to the Milky Way
and the large FOV of the MUSE data mosaic result in a
considerable amount of contamination of nonmember stars in
our Sgr dSph NSC sample. We determine the membership
probability based on the iterative expectation maximization
technique using “clumPy”<fnref n="9" /> (Kimmig et al. 2015). We use our
LOS velocity measurements and the distance of the stars to the
center of the cluster M54 as input values (no metallicity
estimates considered).

This procedure is described in detail by Walker et al. (2009)
and later with some variations by Kimmig et al. (2015). The
technique iteratively estimates the mean velocity and velocity
dispersion of the cluster plus the membership probability of
each star until all the parameters converge. We use a
foreground/background contamination model for the Milky
Way based on the Besançon model of the Galaxy (Robin et al.
2003).

The membership probability obtained with the expectation
maximization technique for each star is represented by color in
the radius versus velocity plot in Figure 3. From our total
sample of 7000 stars we consider those with probability ≥70%

<footnote n="9">⁹ https://github.com/bkimmig/clumpy</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** CMD for the total number of extracted spectra with S/N ≥ 10. Stellar photometry information in F606W (V) and F814W (I) filters (HST/ACS WFC) taken from the M54 catalog (Siegel et al. 2007) of the ACS Survey of Galactic Globular Clusters (Sarajedini et al. 2007). Member stars of the Sgr dSph NSC (probability ≥70%) are color-coded by the S/N logarithm, and nonmembers (probability <70%) are shown in gray.</figcaption></figure>

to be members of the Sgr dSph NSC, leaving a total of 6651 members.

We obtain a maximum likelihood mean LOS velocity for the member stars (dark red points) of $141.34 \pm 0.18 \text{ km s}^{-1}$ and a central velocity dispersion of $16.31 \pm 0.28 \text{ km s}^{-1}$. The quoted numbers are the median $\pm 1\sigma$ values of the posterior distribution. The LOS velocity estimate is in good agreement with previous works, e.g., $V_{LOS} \sim 141 \text{ km s}^{-1}$ by Bellazzini et al. (2008) for the metal-poor and metal-rich populations. From this point, the analysis is carried out on a sample that only includes the members of the Sgr dSph NSC determined by this technique.

In the color–magnitude diagram (CMD) in Figure 4 we show the full sample of extracted single stellar spectra. Member stars of the Sgr dSph NSC (probability $\ge$70%) are color-coded by the S/N logarithm, while nonmembers (probability $<$70%) are shown in gray. The stellar photometry information in F606W (V) and F814W (I) filters is extracted from the M54 HST reference catalog.

### 3.3. Metallicity from Full Spectrum Fitting

The top panel of Figure 5 shows the metallicity histogram for the member stars of the Sgr dSph NSC. This metallicity distribution is highly consistent with the one presented by Mucciarelli et al. (2017) for a total of 76 stars in a radius range of $0.0' < R \le 2.5'$. We cross-matched our sample with those from Carretta et al. (2010a) and Mucciarelli et al. (2017), and we compare the measurements for the stars we have in common. The values are consistent within the errors. We obtain a mean difference of $0.15 \pm 0.03$ dex, when comparing with Carretta et al. (2010a), with a standard deviation of $\sigma = 0.16$. From the comparison with Mucciarelli et al. (2017), we obtain a mean difference of $0.05 \pm 0.02$ dex with a standard deviation of $\sigma = 0.17$. This comparison is included in Appendix B.

In the bottom panel, we present the relation between radius and metallicity for the sample of 6651 stars as derived with ULySS, color-coded by the metallicity uncertainties. Two clear

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Top: metallicity histogram of the member stars of Sgr dSph NSC. Bottom: radius vs. metallicity plot for the same sample. The stars are color-coded by the uncertainties in the [Fe/H] measurements.</figcaption></figure>

distributions in metallicity become apparent: a metal-poor population at [Fe/H] = -1.5 and a metal-rich population around [Fe/H] = -0.3 with a separation at [Fe/H] = -0.8. The median metallicity uncertainty is 0.05 dex.

The top panels in Figure 6 show the CMD color-coded by metallicity (hereafter CMD+metallicity) corresponding to all the Sgr dSph NSC members (and a zoom-in onto the turnoff region; top right panel). We overplotted in blue, green, and orange isochrones from the Dartmouth Stellar Evolution Database (Dotter et al. 2008) corresponding to different stellar populations in the Sgr dSph NSC, using as a reference the work by Siegel et al. (2007). We are able to see stars in different evolutionary stages: turnoff point (TO), red giant branch (RGB), horizontal branch (HB), asymptotic giant branch (AGB), red clump (RC), and even the extreme horizontal branch (EHB) and the blue plume (BP). The CMD shows clear evidence for an old metal-poor population, as well as at least two distinct younger populations with metallicity $\ge$-0.4 that agree well with the 4.5–6 Gyr isochrone highlighted by Siegel et al. (2007).

For our complete sample of member stars we estimate a mean metallicity of [Fe/H] = -1.06, with a standard deviation of $\sigma = 0.64$. For the metal-poor stars ([Fe/H] < -0.8) we obtain a mean metallicity of [Fe/H] = -1.45 and a standard deviation of $\sigma = 0.32$. This is in good agreement with previous works, such as Brown et al. (1999) and Da Costa & Armandroff (1995) with an estimate of [Fe/H] = -1.55, and later by Carretta et al. (2010a), who estimated a metallicity of [Fe/H] = -1.559 $\pm$ 0.021 dex.

In the case of the stars in the metal-rich regime ([Fe/H] > -0.8), we obtain a mean value of [Fe/H] = -0.27 and a
<figure><img src="image_7.png" /><figcaption>Figure 6. Left panels: CMDs for the Sgr dSph NSC members color-coded by metallicity (CMD+metallicity) estimated with ULySS at the top and color-coded by age (CMD+age) at the bottom. Right panels: zoom-in of the TO region. Updated isochrones from the Dartmouth Stellar Evolution Database are overplotted in blue, green, and orange for three stellar populations published by Siegel et al. (2007).</figcaption></figure>

standard deviation of σ = 0.29. In this case, the metallicity measurement differs from previous studies. Carretta et al. (2010a) obtained a value of [Fe/H] = −0.622 ± 0.068 dex, and Mucciarelli et al. (2017) obtained [Fe/H] = −0.52 ± 0.02 dex. A possible explanation for this difference is discussed in Section 4.1.

We repeated the metallicity measurements with a second method, using Ca II triplet lines, and obtained consistent metallicity measurements using the two independent methods (see Appendix C).

3.4. Stellar Age Estimates

Valuable information about the SFH of the Sgr dSph NSC can be obtained from age and metallicity estimates together. For the age estimates of individual stars we use a grid of isochrones, which we compare with the magnitudes and colors of the stars as measured from the HST photometry in V (F606W) and I (F814W) bands with their respective uncertainties (σV, σI), and

metallicities Z, obtained from the MUSE spectra, including their uncertainties σZ. We set the metallicity uncertainty to 0.1 dex when lower than that value. We consider age the only free parameter of the model.

We use Bayes’s theorem to estimate the ages through the posterior probability as

$$
P(\tau|V_{\text{obs}}, I_{\text{obs}}, Z_{\text{obs}}) \propto P(V_{\text{obs}}, I_{\text{obs}}, Z_{\text{obs}}|\tau) \times P(\tau), \quad (1)
$$

where the normalized likelihood function P(Vobs, Iobs, Zobs|τ)
of the observables at a given age is a trivariate Gaussian:

$$
\begin{equation}
\begin{split}
P(V_{\text{obs}}, I_{\text{obs}}, Z_{\text{obs}}|\tau) & \\
& = \frac{1}{\sigma_{V_{\text{obs}}} \sigma_{I_{\text{obs}}} \sigma_{Z_{\text{obs}}} (2\pi)^{3/2}} \exp \frac{-(V_{\text{obs}} - V_0)^2}{2\sigma_{V_{\text{obs}}}} \\
& \quad \times \exp \frac{-(I_{\text{obs}} - I_0)^2}{2\sigma_{I_{\text{obs}}}^2} \exp \frac{-(Z_{\text{obs}} - Z_0)^2}{2\sigma_{Z_{\text{obs}}}^2},
\end{split}
\tag{2}
\end{equation}
$$
where $V_0$, $I_0$, and $Z_0$ denote the real magnitudes and metallicity of the star given its age and the $\sigma$ values represent the independent errors in each measurement. We marginalize over the unknown true magnitudes and metallicity of the star by integrating over the isochrones (see, e.g., Pont & Eyer 2004 for a similar approach).

Due to the lack of extended horizontal branch isochrone models, we have excluded stars with colors between $V-I = -0.2$ and 0.4 and $I$ in the range of 21–17 in our age analysis.

Initially, we used a flat prior in age $P(\tau)$ for all the remaining stars. The age–metallicity relation (AMR) found for the TO region stars turned out to be very well represented by the empirical model published by Layden & Sarajedini (2000). On the other hand, the obtained ages for a fraction of metal-poor RGB stars were lower than expected for such a metallicity regime. This degeneracy could be a consequence of high abundance of elements like oxygen, which has been reported to display a high spread in this metal-poor population (Carretta et al. 2010a). In fact, VandenBerg et al. (2012) showed the impact of different elements abundances on computed evolutionary tracks at similar [Fe/H] values and how this affects the color of the stars and thus their position in the CMD. After this examination of the initial ages, to guard against systematic age uncertainties, we set a weak Gaussian age prior $P(\tau)$ to their model with a standard deviation of 3 Gyr for the age estimates of all extracted member stars from the MUSE cube.

We built the isochrone grid with a resolution of 0.2 Gyr with ages from 0.2 to 15 Gyr from the Dartmouth Stellar Evolution Database (Dotter et al. 2008) assuming the HST/ACS WFC photometric system. The metallicity range is from -2.495 to 0.500 dex with steps of 0.005 dex, considering $[\alpha/Fe] = 0.0$. We correct the isochrone magnitudes by extinction of $A_V = 0.377$ and $A_I = 0.233$ from NED.<fnref n="10" /> These are obtained using the calibration published by Schlafly & Finkbeiner (2011) for the F606W and F814W bands for a reddening of $E(B-V) = 0.15$ estimated in the central 5' radius of M54 (Schlegel et al. 1998) with $R_V = 3.1$. We shifted the grid of isochrones adopting a distance modulus of ($m - M_0$) = 17.27 (Siegel et al. 2007). For comparison, we performed the same analysis considering a reddening of $E(B-V) = 0.16$ and a distance modulus of ($m - M_0$) = 17.13 (Sollima et al. 2010), obtaining consistent age estimates. We take the mode of the likelihood distribution of ages as the best age estimate and the highest density interval that accommodates 68% of the probability as the age uncertainty.

The CMD for the clean stellar sample color-coded by age (hereafter CMD+age) is presented in the bottom left panel of Figure 6, with a zoom-in of the TO region in the bottom right panel of the figure. Three isochrones based on those presented in Siegel et al. (2007) are overplotted; they are also from the Dartmouth Stellar Evolution Database (Dotter et al. 2008).

### 3.5. Populations Split by Age and Metallicity

We use Gaussian Mixture Models (GMMs) in age–metallicity space to separate the stars into distinct multiple populations. To this aim, we used only TO region stars, where we obtain the most reliable age estimates as we describe in

Section 3.4. We tested GMMs with three, four, and five components and found that a GMM with four components gives the best representation of our data set based on the Akaike Information Criterion (Akaike 1974). We optimized the GMM using the affine-invariant Markov chain Monte Carlo (MCMC) algorithm (Goodman & Weare 2010; Foreman-Mackey et al. 2013). We then computed the probabilities of all stars in the observed sample of belonging to any of the GMM defined populations. The fourth Gaussian component of the GMM captures outliers and is not further considered in our analysis.

We estimate the mean and spread in both age and metallicity for each of the subpopulations. We start by taking stars with a high probability (>50%) of belonging to each subpopulation. We then perform a maximum likelihood fit for the ages and metallicities of these stars using a two-dimensional Gaussian and accounting for the observational uncertainties in each stars’ age and metallicity. The result is a mean and *intrinsic* spread in both age and metallicity for each population. We assumed that the [Fe/H] probability distribution functions (pdf’s) are Gaussian and integrated numerically over the non-Gaussian individual stellar age pdf’s. From the results of this procedure we define three different subpopulations as follows:

1. YMR: young metal-rich (red), with mean age of $2.16 \pm 0.03$ Gyr and mean [Fe/H] = $-0.04 \pm 0.01$.

2. IMR: intermediate-age metal-rich (orange), with mean age of $4.28 \pm 0.09$ Gyr and mean [Fe/H] = $-0.29 \pm 0.01$.

3. OMP: old metal-poor (blue), with mean age of $12.16 \pm 0.05$ Gyr and mean [Fe/H] = $-1.41 \pm 0.01$.

We find that the YMR population has an intrinsic $1\sigma$ spread of $0.20 \pm 0.03$ Gyr in age and $0.12 \pm 0.01$ dex in metallicity; the IMR population has an intrinsic spread of $1.16 \pm 0.07$ Gyr in age and $0.16 \pm 0.01$ dex in metallicity; the OMP population has an intrinsic spread of $0.92 \pm 0.04$ Gyr in age and $0.24 \pm 0.01$ dex in metallicity. The uncertainties are represented by the standard deviation in each parameter in the converged part of the MCMC, thus obtaining low errors. The age–metallicity correlation coefficients, $\rho$, for the three populations are $\rho_{YMR} = -0.35$, $\rho_{IMR} = -0.70$, and $\rho_{OMP} = -0.97$.

To check for consistency, we use the same method to estimate the intrinsic spread in metallicity for stars with $I \le 16.5$ mag, where the S/N of the stars is mostly >100. This leaves a total of 22 stars for the YMR (5% of YMR sample) and 208 for the OMP (~9% of OMP sample). For the YMR subpopulation we obtain an intrinsic $1\sigma$ spread of $0.11 \pm 0.06$ dex, and $0.13 \pm 0.03$ dex for the OMP subpopulation. The spread for the YMR subpopulation derived from the brightest stars is consistent with the one measured for the stars of all magnitudes. For the OMP subpopulation, we observe a difference of 0.11 dex. This can be a consequence of the low number of stars in comparison with the entire sample and the effects of the age uncertainties at the RGB region, which can be more uncertain, as we mention in the previous section. Since these values are obtained from a reduced fraction of the entire sample (<10%), we will consider the intrinsic spread measured from the sample at all magnitudes.

The AMR is presented in Figure 7, where we include all stars for which ages were measured and have relative age errors $\le 40\%$. The overplotted crosses represent the intrinsic spreads of the different populations: YMR in red, IMR in orange, and

<footnote n="10">¹⁰ The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Top: age histogram. Bottom: density plot of the age-metallicity relation. The crosses show the intrinsic spread of the different stellar populations in the Sgr dSph NSC: young metal-rich (YMR) in red, intermediate-age metal-rich (IMR) in orange, and old metal-poor (OMP) in blue. The subpopulation's stellar parameters are summarized in Table 1. All the stars included in both panels have age relative errors ≤40%. The vertical dashed lines show where the probability is equal for the two neighboring subpopulations. Note: the inclination of the crosses is given by the correlation coefficient obtained from the Gaussian mixture model.</figcaption></figure>

OMP in blue. The intrinsic age and metallicity spreads, together with the age-metallicity correlation coefficient $\rho$, define the inclination angle of the crosses for each subpopulation. We discuss the possible origin of these populations in Section 5.

For the RC stars in the sample we were able to obtain reliable measurements in metallicity and found that they all fall in the metal-rich regime. The complexity of this stage of stellar evolution prevents us from getting reliable ages for the RC stars from simple isochrone fits. Due to this effect, the RC stars were not considered to belong to any of the subpopulations when we decompose the stars into subpopulations via the maximum likelihood method. For completeness, these stars are shown in Figure 7 and fall on the diagram clustered together as the overdensity at ~1 Gyr in a metallicity range of -0.5 < [Fe/H] < 0.0. To get a better estimate of the likely age distribution of RC stars, we follow Girardi (2016) and compute the RC age distribution for representative SFHs of Sagittarius. We do this for two assumptions of the SFH: (i) the Sgr dSph has a constant SFH, and (ii) the SFH follows the observed SFH from the CMD analysis of de Boer et al. (2015). The de Boer et al. (2015) SFH shows a declining SFR with time but is constructed from the tail of the Sgr dSph, so it is likely missing the most metal-rich and young stars—biasing the star formation to early times. Inclusion of the younger populations in the CMD analysis might produce an SFH more closely represented by a constant SFH. We therefore consider these two cases as plausible bounds on the likely true SFH of Sgr. For the constant SFH case we find that the RC stars are predominantly (~51%) younger than 2.2 Gyr, thus being more likely to belong to the YMR subpopulation than to the IMR (10%). For case (ii), we obtain that anywhere from 10% to 70% of RC stars are younger than 2.2 Gyr, with probabilities of less

than 10% for the IMR. Taking these two cases into consideration, we consider it most likely that the majority of the RC stars are associated with the YMR subpopulation, and we consider them as part of that population in our further analysis.

# 4. Results

## 4.1. Population Analysis

Comparing the metallicity and age information presented in the CMDs in Figure 6, we observe that in the TO region in the CMD+age we distinguish three well-separated subpopulations with different age ranges. The faintest stars have ages over 9 Gyr. This subpopulation belongs to the metal-poor regime in the CMD+metallicity. The middle subgiant branch population has ages of ~3–9 Gyr, while the brightest stars have ages younger than 3 Gyr. These last two populations belong to the metal-rich regime in the CMD+metallicity. From the CMD +metallicity figure, some small difference in metallicity (~0.3 dex) is also visible between these two younger populations.

In the CMD+metallicity the metal-poor RGB is well defined by stars with metallicity below $[Fe/H] = -0.8$, displaying a wide range toward $[Fe/H] = -2.5$, corresponding to an intrinsic iron spread of $\sigma_{[Fe/H]} = 0.24$ dex. In addition, we observe that the stars in the BP region between $18.5 \le I \le 20.5$ and $0.1 \le V - I \le 0.5$ show a wide range in metallicity (from -2.0 to 0.5). We discuss this further in Section 5.3.

From the CMDs and the AMR (Figure 7), we can decipher the SFH of the Sgr dSph NSC. It appears to be extended from 0.5 to 14 Gyr, showing a clear metallicity enrichment toward younger stars.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** For the two panels: red points represent the YMR subpopulation (2.2 Gyr, [Fe/H] = −0.04), orange points are for the IMR (4.3 Gyr, [Fe/H] = −0.29) and blue points are for the OMP (12.2 Gyr, [Fe/H] = −1.41). Gray points show the stars with age relative error greater than 40% or for which age was not estimated. Top panel: CMD. Bottom panel: zoom-in of the top panel showing just the TO region.</figcaption></figure>

We present a CMD comparing the three subpopulations in Figure 8. In both panels the subpopulations are represented as follows: YMR in red, IMR in orange, and OMP in blue. In the top panel we can see how the different populations fall in different positions in the CMD; this is easier to see in the zoom-in of the TO region shown in the bottom panel.

Our findings are in good agreement with the characteristics of three (of four) subpopulations obtained with isochrone fitting by Siegel et al. (2007), and later also by Mucciarelli et al. (2017). From their fitting of the CMD, Siegel et al. (2007) suggest the presence of four populations: (i) the old metal-poor population, with ages of 13 Gyr and [Fe/H] = −1.8; (ii) the intermediate-age metal-rich population, with ages between 4 and 6 Gyr and [Fe/H] = −0.6; (iii) the YMR population with 2.3 Gyr and [Fe/H] = −0.1; and (iv) the very young, very

metal-rich population with [Fe/H] = 0.6 and ages 0.1–0.8 Gyr. We do not detect any stars belonging to their fourth population.

We do find a higher mean metallicity for the metal-rich regime population compared to Carretta et al. (2010a) and Mucciarelli et al. (2017). This may be due to several factors. First, we separate the metal-rich stars into two subpopulations, with the YMR having significantly higher mean metallicity than the IMR; the latter, older stars have metallicity closer to previous measurements. Second, our spatial coverage is confined to much smaller radii in comparison to Carretta et al. (2010a) and Mucciarelli et al. (2017), which cover out to a larger radius. Third, the resolution of our data allows us to extract more stars from the most crowded central regions. Fourth, on average, our metallicity estimates are slightly higher. We obtain a mean $\Delta[\text{Fe/H}]$ of $0.15 \pm 0.03$ dex when comparing with Carretta et al. (2010a), and $0.05 \pm 0.02$ dex when comparing with Mucciarelli et al. (2017). The differences in the measurements might be caused by the method used and/or spectral resolution (see Appendix B).

In Table 1 we include a summary of all the estimated parameters for the subpopulations.

## 4.2. Subpopulation Spatial Distributions

We now analyze the spatial distribution of the three subpopulations in the Sgr dSph NSC defined in the previous section. In Figure 9 we present the cumulative radial distribution for the three subpopulations: YMR (red), IMR (orange), and OMP (blue). We consider all the member stars extracted with $S/N \ge 10$. In order to fairly compare the three subpopulations’ distribution in terms of completeness, we also constrain the sample to magnitudes $I \le 20.5$ mag.

In Figure 9 we observe that the YMR subpopulation is the most centrally concentrated of the three in this NSC. The second-highest central density is shown by the OMP subpopulation, which is the dominant one in stellar number with over 2000 stars. Finally, the lowest central concentration is shown by the IMR subpopulation. Comparing the distribution of these subpopulations with a uniform distribution (magenta dashed line), we observe that the stars in the IMR subpopulation are the least centrally concentrated in the most central $40''$. However, they are still not uniformly distributed; as we observed in Figure 9, the distribution shows that the stars in the IMR subpopulation are still significantly centrally peaked. This difference in the spatial distribution between the subpopulations suggests different origins. We discuss this in detail in Section 5.

It is important to consider that, due to extreme central crowding and the limitations of our spatial resolution, we are not able to extract all the stars that actually reside in this NSC. Improved spatial resolution is needed to more accurately derive the spatial distribution of the different stellar populations.

## 4.3. 2D Morphology of Subpopulations

Using the coordinate information of the stars and the data FOV, we fit a two-dimensional Plummer profile (Plummer 1911) to each of the subpopulations, optimizing the best profile using an MCMC algorithm (Goodman & Weare 2010; Foreman-Mackey et al. 2013). We set the centroid, ellipticity, position angle (PA), and half-light radius as free parameters, thus obtaining the best estimates for each stellar subpopulation.
<table><caption>**Table 1**
Summary of the Stellar Subpopulations in M54</caption><thead><tr><th>Subpopulations</th><th>YMR</th><th>IMR</th><th>OMP</th></tr></thead><tbody><tr><td>[Fe/H]</td><td>-0.04 ± 0.01</td><td>-0.29 ± 0.01</td><td>-1.41 ± 0.01</td></tr><tr><td>σ<sub>[Fe/H]</sub></td><td>0.12 ± 0.01</td><td>0.16 ± 0.01</td><td>0.24 ± 0.01</td></tr><tr><td>Age (Gyr)</td><td>2.16 ± 0.03</td><td>4.28 ± 0.09</td><td>12.16 ± 0.05</td></tr><tr><td>σ<sub>Age</sub> (Gyr)</td><td>0.20 ± 0.03</td><td>1.16 ± 0.07</td><td>0.92 ± 0.04</td></tr><tr><td>Correlation factor (ρ)</td><td>-0.35</td><td>-0.70</td><td>-0.97</td></tr><tr><td>Median V<sub>LOS</sub> (km s<sup>-1</sup>)</td><td>141.92 ± 0.54</td><td>142.61 ± 0.59</td><td>141.22 ± 0.26</td></tr><tr><td>R.A. centroid (deg)</td><td>283.76353 ± 0.00174</td><td><sup>a</sup></td><td>283.76299 ± 0.00094</td></tr><tr><td>Decl. centroid (deg)</td><td>-30.477006 ± 0.000986</td><td><sup>a</sup></td><td>-30.477054 ± 0.000680</td></tr><tr><td>Half-light radius (arcmin)</td><td>1.47 ± 0.20</td><td><sup>a</sup></td><td>1.90 ± 0.12</td></tr><tr><td>Ellipticity</td><td>0.31 ± 0.10</td><td><sup>a</sup></td><td>0.16 ± 0.06</td></tr><tr><td>Position angle (degrees)</td><td>4.56 ± 10.84</td><td><sup>a</sup></td><td>16.83 ± -12.81</td></tr><tr><td>Number of stars</td><td>440</td><td>536</td><td>2550</td></tr></tbody></table>

**Note.**

<sup>a</sup> Measurements do not converge. σ<sub>[Fe/H]</sub> and σ<sub>Age</sub> correspond to the [Fe/H] and age intrinsic spreads, respectively. V<sub>LOS</sub>: line-of-sight velocity. The correlation factor (ρ) indicates how age and metallicity are correlated and is estimated using Gaussian Mixture Models.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Cumulative radial distribution for the three subpopulations for the member stars extracted with S/N ≥ 10 and magnitudes I ≤ 20.5 mag. The YMR subpopulation is shown in red (2.2 Gyr, [Fe/H = −0.04], IMR in orange (4.3 Gyr, [Fe/H] = −0.29), and OMP in blue (12.2 Gyr, [Fe/H] = −1.41). The magenta dashed line describes a uniform stellar distribution. The vertical black dashed line shows the half-light radius of M54 (R<sub>HL</sub> = 0.82; Harris 1996, 2010 edition).</figcaption></figure>

The estimated parameters are summarized in Table 1. The likelihoods for the YMR and OMP populations converge and give good estimates for the centroid, ellipticity, and PA. However, for the IMR it does not converge owing to the less centrally concentrated distribution that this subpopulation shows (see Section 4.2).

## 5. Discussion

We show the systemic velocities of the three subpopulations in Table 1. The difference of ~1σ to 2σ between them could be a consequence of the number of stars in each sample, which is considerably higher for the OMP, with ~5 times more stars than the YMR and IMR subpopulations. This results in lower errors for the OMP. In spite of the slight difference, the velocities are consistent. The nearly identical systemic velocities derived from a high number of stars strongly support the idea that all three subpopulations are spatially coincident, as also suggested in

previous studies (e.g., Da Costa & Armandroff 1995; Ibata et al. 1997; Monaco et al. 2005a; Bellazzini et al. 2008). We also add another piece of evidence to this with their respective centroids (Table 1), showing that the subpopulations OMP and YMR are actually spatially coincident with the same centroid. This strongly argues against a chance alignment in projection between the metal-poor and metal-rich populations (Siegel et al. 2011).

### 5.1. A Possible Merger Remnant as the Seed of the Sgr dSph NSC

We find that the stars in the metal-poor regime show a large spread in age, σ<sub>age</sub> = 0.9 Gyr, and iron content, σ<sub>[Fe/H]</sub> = 0.24, a higher spread in comparison with the literature (σ<sub>[Fe/H]</sub> = 0.186; Carretta et al. 2010b; Willman & Strader 2012). Large iron spreads in GCs can be explained by a scenario in which two GCs merge (Bekki & Tsujimoto 2016; Gavagnin et al. 2016; Khoperskov et al. 2018, and references therein).

Gavagnin et al. (2016), using N-body simulations, studied the structural and kinematic signatures in the remnants of GC mergers with a sample of progenitors of different densities and masses. They pointed out that if two GCs were formed in the same dwarf galaxy (or molecular cloud) and display a low relative velocity, they could merge, making the merging scenario more likely in dwarf galaxies. They actually mention the Sgr dSph as a good candidate for this scenario to happen. Bellazzini et al. (2008) found that the velocity dispersion profile is close to flat, with values of ~10 km s<sup>-1</sup> in the innermost 80'.

Bekki & Tsujimoto (2016) used numerical simulations to suggest that GC mergers explain the existence of complexity in high-mass GCs, with multiple stellar populations and a large spread in metallicity, what the authors called “anomalous” GCs. Performing several tests, Bekki & Tsujimoto (2016) found that the merger between GCs with masses greater than 3 × 10<sup>5</sup> M<sub>☉</sub> will occur inevitably in the host dwarf galaxy, with a mass in the range of M<sub>dh</sub> = 3 × 10<sup>9</sup>–3 × 10<sup>10</sup> M<sub>☉</sub>. This is due to stronger dynamical friction effects on the GCs compared to the field stars in those shallow, low-dispersion galaxies. The authors mention that if the clusters are massive enough, the timescale for the merger is just a few gigayears, and it occurs before the total disruption of the dwarf galaxy due to
the Galactic tidal field. They find that with a merger of at least two GCs it would be possible to observe GCs with an internal metallicity spread. In the case of the Sgr dSph, the progenitor dark halo mass before infall to the Milky Way has been estimated with a lower limit of $6 \times 10^{10} M_{\odot}$ by Gibbons et al. (2017), in good agreement with Mucciarelli et al. (2017). The mass estimated for M54 is $(1.41 \pm 0.02) \times 10^6 M_{\odot}$ (Baumgardt & Hilker 2018), consistent with mergers of GCs of $10^5 M_{\odot}$.

The merger scenario might explain the high mass and the large iron and age spreads in the OMP subpopulation. The mass of the Sgr dSph progenitor would be conducive for the GCs to merge. From simulations, the infall of a massive GC to the center of the Sgr dSph due to dynamical friction effects has been estimated to have occurred ~5–9 Gyr ago (Bellazzini et al. 2006, 2008). The infall time is at most 3 Gyr (Bellazzini et al. 2008) for a cluster with M54’s mass. The spread in age and metallicity of the oldest, metal-poor stars is consistent with the scenario where multiple GCs were driven to the center of the Sgr dSph and merged, thus building up the metal-poor component of the NSC.

## 5.2. In Situ Formation in the Sgr dSph NSC: YMR Subpopulation

From the metallicity and age characterization of our stellar sample we distinguish two subpopulations in the metal-rich regime. We are able to separate them using the AMR and CMD position of the stars. We observe that the youngest population —YMR—has the highest central density in the Sgr dSph NSC. The high metallicity, young age, and spatial distribution suggest that the formation of the YMR subpopulation formed in situ starting ~3 Gyr ago.

Additional evidence for in situ formation of the YMR population is that it appears more flattened than the OMP population. We measure an ellipticity of $0.31 \pm 0.10$ for the YMR component, while the OMP population has an ellipticity of $0.16 \pm 0.06$ (Table 1). Centrally concentrated and sometimes flattened young subcomponents are commonly found in other NSCs in both early- and late-type galaxies (Seth et al. 2006; Carson et al. 2015; Feldmeier-Krause et al. 2015; Nguyen et al. 2017); given the young age of these populations, in situ formation is strongly favored, with the flattening suggestive of formation in a star-forming gas disk. We note that the difference in ellipticity between the YMR and OMP subpopulations is small, given the measurement uncertainties. However, this difference becomes meaningful when we combine the ellipticity with the kinematics. In a second paper we will present the kinematic characterization of the subpopulations, adding a new diagnostic to address this issue.

These young, flattened populations can survive despite being embedded in an older, hotter population. Mastrobuono-Battisti & Perets (2013, 2016) modeled the evolution of stellar disks in dense stellar clusters using N-body simulations, finding that even after several gigayears (around the Milky Way GC’s age) the stars are still not fully mixed. Thus, both generations of stars show different distributions. They observed that the second population is concentrated at the center after 12 Gyr with no relaxed spherical shape. The amount of rotation is also evidence for this model, where the second generation of stars is found to rotate faster than the first one. From our kinematic analysis, we found that the YMR subpopulation rotates faster than the OMP. We will present the evidence on this matter in our upcoming paper on this NSC’s kinematics.

To form the youngest metal-rich subpopulation, enriched gas is needed; however, no neutral gas has been detected in this galaxy. Using simulations, Tepper-García & Bland-Hawthorn (2018) found that the Sgr dSph might have lost all its gas after the second encounter with the Galactic disk. They found that two pericentric passages happened ~2.8 and ~1.3 Gyr ago. The first produced a gas loss of 30%–50%, while the second would be responsible for the stripping of all the residual gas. The timescales of both encounters are consistent with the latest bursts of star formation seen in Figure 7 and were suggested based solely on the photometry by Siegel et al. (2007). This suggests that the YMR population may have formed by gas being driven into the nucleus during encounters with the Milky Way, supporting the in situ formation of this population.

The in situ formation of a new generation of stars could alternatively be explained for a massive NSC with a deep enough potential well to collect the enriched gas ejected from metal-rich and high-mass stars at larger radii that cools and sinks to the center (Bailey 1980). This seems plausible given the characteristics of the stars in the IMR subpopulation and the relatively small fraction of YMR stars in comparison with the OMP.

## 5.3. Sgr dSph Field Stars in Its NSC: IMR Subpopulation

A different situation is observed for the IMR population in comparison to the YMR. These stars seem to be less centrally concentrated in the field, with a wide spread in ages from 3 to 8 Gyr. These findings suggest that the star formation in the Sgr dSph during this period was not particularly concentrated in the nucleus. The metallicity range of these stars is consistent with those found for stars close to the center of the galaxy (Hasselquist et al. 2017). In addition, comparing the IMR sample of our AMR (Figure 7) with the one presented by de Boer et al. (2015) for the Sagittarius bright stream, we see a good agreement in both age (4–14 Gyr) and metallicity ($-2.5 < [Fe/H] < -0.1$), thus supporting our good recovery of the underlying chemical evolution of the host galaxy.

The spatial distribution of these two subpopulations—YMR and IMR—has also been reported by Mucciarelli et al. (2017). Using a sample of 109 stars with $[Fe/H] \ge -1.0$ and a magnitude limit of $I_{mag} = 18$, the authors observed a shift in the metallicity peak in the metal-rich population at different projected distances from its center. At $0.1' < R < 2.5'$ the peak is at $[Fe/H] = -0.38$, changing to $[Fe/H] = -0.45$ at $2.5' < R < 5.0'$, noticing a metallicity gradient for this population. In addition, the authors present the cumulative radial distribution of the two young populations, where they found that the youngest subpopulation is more centrally concentrated than the intermediate-age population. With this finding they suggest that the youngest population is the dominant population in the Sgr dSph NSC, with the intermediate-age one becoming more important at larger radii. Since our observations reach out to ~2.0' from the center of this NSC, we are not able to see a change in the peak at larger radii. However, we also see this behavior between the two youngest metal-rich subpopulations.

Our spectroscopic sample also includes stars in the BP region, which could be populated by YMR stars of $[Fe/H] = 0.6$ and $0.1-0.8$ Gyr (Siegel et al. 2007) or blue straggler stars (BSSs). However, we measure a wide range of metallicities for these stars, suggesting that they are BSSs instead of young stars, which would be expected to display a more homogeneous
metallicity, e.g., as observed in the YMR subpopulation. Mucciarelli et al. (2017) found that the cumulative radial distributions of the BP and the intermediate-age metal-rich stars were not distinguishable. Since this is a less dense environment in comparison to the center of the NSC, these BSSs could be the product of mass transfer between binaries as has been found in other dSphs (e.g., Mapelli et al. 2007, 2009; Momany et al. 2007).

As discussed in this section, the stars from the IMR subpopulation seem to have properties consistent with those at larger radii in the galaxy. This population shows a central concentration, shallower than the other two. With the current information we are not able to tell if the stars in the inner regions are actually dynamically bound to the NSC (with the YMR and OMP) or if they are part of the main body of the host galaxy.

## 5.4. The Formation History of the Sgr dSph NSC

From this work and previous studies of the populations in the Sgr dSph NSC, we can put together the story of this NSC. It starts with two or more massive GCs that eventually merge at the center of the Sgr dSph, forming a massive nucleus consisting of old and metal-poor stars with a large metallicity spread. The two encounters of the Sgr dSph with the Galactic disk occurred ~2.8 and ~1.3 Gyr ago (Tepper-García & Bland-Hawthorn 2018) and could have triggered two new episodes of star formation before the total stripping of the gas, creating in the first the YMR subpopulation. This results in a complex multipopulation NSC.

This NSC is on its way to becoming a stripped nucleus considering the ongoing strong interaction between its host—the Sgr dSph—and the Milky Way. In fact, Bellazzini et al. (2008) suggested that this nucleus will probably end up as a compact remnant with two populations: old metal-poor and young metal-rich, with no signatures of the progenitor galaxy. This puts the Sgr dSph NSC in close company with the most massive cluster in the Milky Way, ω Cen, which has long been considered a potential stripped nucleus (e.g., Lee et al. 1999; Carretta et al. 2010b). ω Cen presents a centrally concentrated disk-like component (van de Ven et al. 2006, see their Figures 19 and 20), very similar to the YMR subpopulation we detect in the Sgr dSph NSC. The age spread in ω Cen (at least 2 Gyr; Hilker et al. 2004; Villanova et al. 2014), similar to the spread we see in the OMP, suggests the merging of GCs early on in the nucleus of a progenitor galaxy (e.g., Bekki & Tsujimoto 2016). Unlike ω Cen, the stripping around the Sgr dSph NSC is ongoing, and thus we have the opportunity to understand the role mergers and stripping have played in creating the cluster we see today.

Given all the evidence we have presented in this paper, our suggestion to the community is to revert to the original naming and use “M54” in the same manner as it was given to the object upon its discovery by Charles Messier in 1778. The evidence consistently suggests that M54 is not a normal metal-poor GC but a complex NSC.

## 6. Conclusion

In this work we present a rich sample of ~6600 stellar spectra extracted from a mosaic of 16 pointings of MUSE data on M54, the NSC of the Sgr dSph, a dwarf galaxy currently

being disrupted by the Milky Way. Through LOS velocity, metallicity ([Fe/H]), and age measurements we have characterized M54's stellar populations. We were able to detect at least three subpopulations with the same systemic velocity, differentiated by age and metallicity, where two of them have the same centroid.

The subpopulations we find are as follows: (i) YMR: young metal-rich, with ages 2.2 Gyr and average metallicity [Fe/H] = -0.04; (ii) IMR: intermediate-age metal-rich, with ages of 4.3 Gyr and metallicity [Fe/H] = -0.29; and (iii) OMP: old metal-poor, with ages 12.2 Gyr and metallicity [Fe/H] = -1.41.

The existence of these three distinct subpopulations with the displayed differences in age and metallicity suggests the following conclusions:

1. The stars in the OMP population have ages and metallicity consistent with it being assembled by two or more star clusters in spiraling to the nucleus via dynamical friction.

2. The YMR population is both more flattened and more centrally concentrated than the other two populations. These features suggest in situ formation from enriched gas retained in the deep potential well of M54. This young, centrally concentrated component is similar to features observed in other NSCs. We estimate that the YMR subpopulation formation episode started around 3 Gyr ago, consistent with the time of the first big encounter between the Sgr dSph and the Milky Way, suggesting that gas was channeled into the nucleus during this encounter. The youngest stars in the YMR population, <3 Gyr, might be related to when the Sgr dSph lost its gas during its ongoing interaction with the Milky Way. Alternatively, the YMR subpopulation could have formed from gas ejected from high-mass and metal-rich stars in the IMR subpopulation retained in the deep potential of the massive OMP subpopulation.

3. Our metallicity measurements for the IMR subpopulation are consistent with those for the field star population of the galaxy, including regions close to the center. This subpopulation shows the lowest central concentration in M54 but is still significantly centrally peaked. Additional information is needed to determine whether these stars are actually dynamically bound to the NSC.

M54 is a unique test case. In this complex nucleus we find evidence for two processes that build up the NSC: (i) infall of two or more GCs, which merge to create a single high-mass cluster with a large metallicity spread, and (ii) in situ star formation from enriched gas in the nucleus. In this case, the first scenario could be the key for the second to occur. This detailed formation history of the Sgr dSph NSC helps us understand the processes of NSC formation and the role of galaxy-galaxy interaction in this formation.

In a second paper we will present a kinematic analysis for the sample. We will include a kinematic analysis for the three age-metallicity subpopulations found in this work, giving additional insight on M54's formation.

We thank the referee, Michele Bellazzini, for the constructive comments and suggestions that helped to improve this work. A.C.S acknowledges financial support from NSF grant AST-1350389. S.K. gratefully acknowledges funding from a European Research Council consolidator grant (ERC-CoG-
646928- Multi-Pop). R.L. acknowledges funding from a Natural Sciences and Engineering Research Council of Canada PDF award, and R.L. and N.N. acknowledge support by Sonderforschungsbereich SFB 881 "The Milky Way System" (subprojects A7, A8, and B8) of the German Research Foundation (DFG). M.A.-C., R.L., and P.B. acknowledge support from DAAD PPP project No. 57316058 "Finding and exploiting accreted star clusters in the Milky Way." G.v.d.V. and L.L.W. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 724857 (Consolidator Grant ArcheoDyn). This research has made use of NASA's Astrophysics Data System. We thank Brian Kimmig for kindly providing his code "clumPy" (see footnote 9) for stellar cluster membership estimates. We thank Alessandra Mastrobuono-Battisti and Iskren Georgiev for their helpful comments and discussions on this work.

# Appendix A

## Comparing Velocity and Metallicity Measurements for Multiple Stellar Spectra

In Section 2.2 we describe the stellar spectra extraction from the MUSE data set. Due to overlaps between the neighboring pointings, we have from one up to four spectra for an individual star. To compare how the velocity, metallicity, and their respective uncertainties are correlated in this appendix, we

include the comparison between the velocity estimates for the complete extracted stellar sample (see Section 3.1). The top left panel in Figure 10 shows velocity versus velocity measurements of the repeated stars. We show the same for their corresponding errors in the middle left panel. The bottom left panel shows the histogram of the difference in the measurements. The estimated sigma is $\sigma = 16.25 \text{ km s}^{-1}$.

The right panels of Figure 10 include the same comparison for the case of the metallicity measurements described in Section 3.3 for the same stellar sample. The top right panel of Figure 10 shows metallicity versus metallicity measurements of the repeated stars. The correlation between their respective errors is shown in the middle right panel. The bottom left panel shows the histogram of the difference in the measurements, with a $\sigma$ of 0.51 dex.

Since brighter stars have higher S/N, to show the variation of the spread in dependency with magnitude, in Table 2 we include the standard deviation of the [Fe/H] and LOS velocity measurements at different magnitude cuts.

The difference in the values for both parameters can be due to the different observing conditions during the acquisition of the 16 pointings (between 2015 June and July). This produces variations in the quality of the spectra. In spite of this, we extracted good-quality spectra with a reasonable S/N threshold of 10. The dispersion we observe between repeated stars is acceptable and does not affect our main results.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Top left panel: velocity measurements' correlation of the repeated stars. Middle left panel: correlation between the uncertainties. Bottom left panel: histogram of the difference in the measurements, with $\sigma = 16.25$ km s⁻¹. Top right panel: metallicity measurements' correlation of the repeated stars. Their respective errors' correlation is shown in the middle right panel. Bottom left panel: histogram of the difference in the measurements, with $\sigma = 0.51$ dex.</figcaption></figure>
<table><caption>**Table 2**
Line-of-sight Velocity and [Fe/H] Spread for Repeated Stars</caption><thead><tr><th>Imag Range</th><th>σ<sub>Vr</sub> (km s<sup>-1</sup>)</th><th>σ<sub>[Fe/H]</sub></th><th>Number of Stars</th></tr></thead><tbody><tr><td>[13–14]</td><td>5.84</td><td>0.19</td><td>55</td></tr><tr><td>[14–15]</td><td>12.75</td><td>0.40</td><td>135</td></tr><tr><td>[15–16]</td><td>13.12</td><td>0.51</td><td>235</td></tr><tr><td>[16–17]</td><td>16.83</td><td>0.54</td><td>425</td></tr><tr><td>[17–18]</td><td>17.78</td><td>0.54</td><td>430</td></tr><tr><td>[18–19]</td><td>12.30</td><td>0.47</td><td>324</td></tr><tr><td>[19–20]</td><td>16.87</td><td>0.50</td><td>236</td></tr><tr><td>[20–21]</td><td>29.51</td><td>0.58</td><td>68</td></tr><tr><td>[21–22]</td><td>30.75</td><td>0.42</td><td>6</td></tr><tr><td>[13–22]</td><td>16.25</td><td>0.51</td><td>1 920</td></tr></tbody></table>

Note. σ corresponds to the spread given by the standard deviation of the measurements from the repeated stars at different magnitude cuts.

## Appendix B
Metallicity Comparison with Previous Works

The FOV of our MUSE data partially overlaps with those in Carretta et al. (2010a) and Mucciarelli et al. (2017). We cross-matched the samples using coordinate information for the stars with

$I \le 17.75$ mag, since the S/N of the spectra is larger and warrants reliable metallicity measurements. In the top panel of Figure 11 we present the difference in [Fe/H] ($\Delta$[Fe/H]) of the 51 stars that are in common between this work and Mucciarelli et al. (2017) (red circles), and the 36 with Carretta et al. (2010a) (blue circles). The dashed black line shows a $\Delta$[Fe/H] = 0. For $\Delta$[Fe/H] between this work and Mucciarelli et al. (2017), we obtain a mean value of $0.05 \pm 0.02$ dex, with a standard deviation of $\sigma = 0.17$. For $\Delta$[Fe/H] between this work and Carretta et al. (2010a), we have a mean value of $0.15 \pm 0.03$ dex, with a standard deviation of $\sigma = 0.16$. With the respective colors, we include in the figure the $\Delta$[Fe/H] mean values as solid lines; the dashed lines show the area confined by the error of the mean. In this analysis we observe a good agreement between our metallicity measurements and those of previous studies, which use different methods and, in the case of Carretta et al. (2010a), high-resolution spectra.

In the bottom panel of Figure 11 we show the relation between the metallicity values from this work and those of Mucciarelli et al. (2017; 51 stars; red circles) and Carretta et al. (2010a; 36 stars; blue circles). In addition, we include the 13 stars in common between Carretta et al. (2010a) and Mucciarelli et al. (2017) as gray circles. The cross-match between these 13 stars and our sample gives eight stars in

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

**Figure 11.** Top: difference in metallicity [Fe/H] between this and other works for stars with $13 \le I \le 17.75$ mag. The 51 red circles show the $\Delta$[Fe/H] between this work and Mucciarelli et al. (2017), and the 36 blue circles are for Carretta et al. (2010a). The mean $\Delta$[Fe/H] between this work and Mucciarelli et al. (2017) is $0.05 \pm 0.02$ dex, with $\sigma = 0.17$. For $\Delta$[Fe/H] between this work and Carretta et al. (2010a), the mean is $0.15 \pm 0.03$ dex, with $\sigma = 0.16$. The solid red and blue lines show the mean value of the respective $\Delta$[Fe/H] estimates, while the same color dashed lines show the area confined by the mean error. The dashed black line shows a $\Delta$[Fe/H] = 0. Bottom: relation between [Fe/H] estimates from this work and those of other works: red circles for Mucciarelli et al. (2017; 51 stars), and blue circles for Carretta et al. (2010a; 36 stars). The gray circles represent the relation of 13 common stars between Carretta et al. (2010a) and Mucciarelli et al. (2017). The eight common stars between the three samples are represented as green circles when comparing this work with Mucciarelli et al. (2017) and green squares when comparing with Carretta et al. (2010a). The dashed black line shows a one-to-one relation.
common. These stars are included in green circles when comparing this work with Mucciarelli et al. (2017) and green squares when comparing with Carretta et al. (2010a).

## Appendix C
Ca II Triplet Metallicity

For the metallicity estimates, we use two independent methods to test the robustness of our results. The second method we use to estimate the metallicity is based on the calcium triplet lines (Ca II triplet) at 8498, 8542, and 8662 Å (see, e.g., Starkenburg et al. 2010; Carrera et al. 2013). We use the calibration presented by Carrera et al. (2013), obtained using the Ca II triplet and calibration measurements obtained from hundreds of RGB star observations.

These nonlinear calibrations are specifically tailored toward the metal-poor regime, reaching reliable estimates of metallicity for stars between -4.0 < [Fe/H] < 0.5. Ca II triplet lines are strong on the RGB, but their strength depends on the joint effects from metallicity, temperature, and surface gravity. Therefore, to obtain just the effect due to the metallicity, the authors removed these last two contributions using luminosity indicators in order to obtain the following relation (Equation (2) in Carrera et al. 2013):

$$ [\text{Fe/H}] = a + b \times \text{Mag} + c \times \Sigma\text{Ca} + d \times \Sigma\text{Ca}^{1.5} + e \times \Sigma\text{Ca} \times \text{Mag.} \quad (3) $$

We use Mag—the luminosity indicator—as V-V_HB, where V_HB = 18.16 (Harris 1996) is the average magnitude of the HB. The coefficient values according to Carrera et al. (2013) are as follows: a = -3.45, b = 0.11, c = 0.44, d = -0.65, e = 0.03. ΣCa corresponds to the sum of the equivalent width (ΣEW) of the three Ca II lines. For additional details see Carrera et al. (2013).

We applied the Doppler correction to the spectra using the LOS velocity values of each star and performed a continuum normalization. We fit a Gauss and Voigt profile to estimate the EW of the Ca II triplet lines. We use the spectral ranges for these lines defined by Armandroff & Zinn (1988). We find that Voigt profiles provided better fits than Gaussian profiles for high S/N; therefore, we have used Voigt profile fits for estimating the metallicity. In Figure 12 we show fitting examples for two stars with different S/N values (S/N = 10 and S/N = 95). The best-fit Voigt and Gaussian profiles are overplotted as dashed blue and solid red lines, respectively.

Using the ΣEW values from the Voigt fitting in the calibration in Equation (3), we obtain metallicity values from -2.9 < [Fe/H] < 0.5. Since this calibration is limited to RGB stars, for this procedure we consider a cut in V-I = 0.7, selecting 4051 stars for the metallicity estimates in the sample of members of the Sgr dSph NSC. In the top panel of Figure 13 we show the ΣEW versus V-V_HB; overplotted in red lines are the metallicity values for this calibration. The middle panel shows the metallicity histogram of the metallicity values obtained with this method (red). For comparison, we have overplotted the metallicity histogram of the estimates obtained performing full spectral fitting using ULySS (Koleva et al. 2009; blue) for the same sample of RGB stars (see Section 3.3). The vertical black line on the left shows the average metallicity

<figure><img src="image_13.png" /><figcaption>Figure 12. Fitting examples of Voigt (dashed blue line) and Gaussian (solid red line) profiles to the Ca II triplet lines for two spectra. The dashed gray line shows the observed spectrum. The black solid line is the selected region to evaluate the continuum level. Red crosses mark show the discarded pixels from the continuum region. Top: example of high-S/N spectrum (S/N = 95). Bottom: example of low-S/N spectrum (S/N = 10).</figcaption></figure>

for the metal-poor population ([Fe/H] = -1.559 ± 0.021 dex), estimated by Carretta et al. (2010a) using neutral Fe lines of 76 stars. The vertical black line to the right shows the average metallicity for 25 metal-rich stars ([Fe/H] = -0.622 ± 0.068 dex) estimated by the same authors. We see consistent results between our metallicity measurements using full spectral fitting with ULySS (Koleva et al. 2009) and the measurement by Carretta et al. (2010a) for the metal-poor stars. For the metal-rich population, we see a difference in the average metallicity, where our value is more metal-rich ([Fe/H] ~ -0.3 dex). We give a possible explanation for this difference in Section 4.1. The bottom panel shows the correlation of the metallicity values between the two methods, where we see a bigger difference in the high-metallicity regime. The estimates are plotted in blue for the spectra with S/N > 60 and in gray for S/N < 60.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Top panel: equivalent width sum of the three Ca II triplet lines vs. $V - V_{HB}$. Overplotted in red lines are the metallicity values for the Carrera calibration using the spectral ranges of the Ca II triplet lines defined by Armandroff & Zinn (1988). Middle panel: metallicity histogram of the values obtained with the Ca II triplet line method using the calibration published by Carrera et al. (2013) (red) and those obtained performing full spectral fitting using ULySS (Koleva et al. 2009) for the same stellar sample (blue). There are fewer stars in the first method since it is restricted to RGB stars (~4000). From left to right the vertical black lines show the average metallicity for the metal-poor and metal-rich stars estimated by Carretta et al. (2010a). Bottom panel: correlation of the metallicity values obtained with the two methods, for spectra with S/N > 60 in blue and with S/N < 60 in gray.</figcaption></figure>

## ORCID iDs

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

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

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

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

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

## References

Afanasiev, A. V., Chilingarian, I. V., Mieske, S., et al. 2018, MNRAS, 477, 4856

Ahn, C. P., Seth, A. C., Cappellari, M., et al. 2018, ApJ, 858, 102

Ahn, C. P., Seth, A. C., den Brok, M., et al. 2017, ApJ, 839, 72

Akaike, H. 1974, ITAC, 19, 716

Armandroff, T. E., & Zinn, R. 1988, AJ, 96, 92

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

Bailey, M. E. 1980, MNRAS, 191, 195

Bailin, J. 2018, ApJ, 863, 99

Baumgardt, H. 2017, MNRAS, 464, 2174

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

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

Bellazzini, M., Correnti, M., Ferraro, F. R., Monaco, L., & Montegriffo, P. 2006, A&A, 446, L1

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

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

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

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

Brown, J. A., Wallerstein, G., & Gonzalez, G. 1999, AJ, 118, 1245

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

Carrera, R., Pancino, E., Gallart, C., & del Pino, A. 2013, MNRAS, 434, 1681

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

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

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

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

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

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

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

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

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

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

Ferraro, F. R., Massari, D., Dalessandro, E., et al. 2016, ApJ, 828, 75

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

Fuentes-Carrera, I., Jablonka, P., Sarajedini, A., et al. 2008, A&A, 483, 769

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

Genzel, R., Eisenhauer, F., & Gillessen, S. 2010, RvMP, 82, 3121

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

Georgiev, I. Y., Hilker, M., Puzia, T. H., Goudfrooij, P., & Baumgardt, H. 2009, MNRAS, 396, 1075

Georgiev, I. Y., Neumayer, N., Lützgendorf, N., Böker, T., & Leigh, N. 2016, MNRAS, 457, 2122

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

Girardi, L. 2016, ARA&A, 54, 95

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

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

Hasselquist, S., Shetrone, M., Smith, V., et al. 2017, ApJ, 845, 162

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

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

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

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

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

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

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

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

Kimmig, B., Seth, A., Ivans, I. I., et al. 2015, AJ, 149, 53
Koleva, M., Prugniel, P., Bouchard, A., & Wu, Y. 2009, A&A, 501, 1269

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

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

Majewski, S. R., Skrutskie, M. F., Weinberg, M. D., & Ostheimer, J. C. 2003, ApJ, 599, 1082

Mapelli, M., Ripamonti, E., Battaglia, G., et al. 2009, MNRAS, 396, 1771

Mapelli, M., Ripamonti, E., Tolstoy, E., et al. 2007, MNRAS, 380, 1127

Massari, D., Mucciarelli, A., Ferraro, F. R., et al. 2014, ApJ, 795, 22

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

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

Meylan, G., Sarajedini, A., Jablonka, P., et al. 2001, AJ, 122, 830

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

Momany, Y., Held, E. V., Saviane, I., et al. 2007, A&A, 468, 973

Monaco, L., Bellazzini, M., Bonifacio, P., et al. 2005a, A&A, 441, 141

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

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

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

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

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

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

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

Ordenes-Briceño, Y., Puzia, T. H., Eigenthaler, P., et al. 2018, ApJ, 860, 4

Origlia, L., Mucciarelli, A., Fiorentino, G., et al. 2019, ApJ, 871, 114

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

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

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

Pont, F., & Eyer, L. 2004, MNRAS, 351, 487

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

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

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

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Schödel, R., Feldmeier, A., Kunneriath, D., et al. 2014a, A&A, 566, A47

Schödel, R., Feldmeier, A., Neumayer, N., Meyer, L., & Yelda, S. 2014b, CQGra, 31, 244007

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

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

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

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

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

Sohn, S. T., van der Marel, R. P., Carlin, J. L., et al. 2015, ApJ, 803, 56

Sollima, A., Cacciari, C., Bellazzini, M., & Colucci, S. 2010, MNRAS, 406, 329

Starkenburg, E., Hill, V., Tolstoy, E., et al. 2010, A&A, 513, A34

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

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

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

VandenBerg, D. A., Bergbusch, P. A., Dotter, A., et al. 2012, ApJ, 755, 15

Vasiliev, E. 2019, MNRAS, 484, 2832

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

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

Walker, M. G., Mateo, M., Olszewski, E. W., Sen, B., & Woodroofe, M. 2009, AJ, 137, 3109

Weilbacher, P. M., Streicher, O., Urrutia, T., et al. 2014, in ASP Conf. Ser. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay (San Francisco, CA: ASP), 451

Willman, B., & Strader, J. 2012, AJ, 144, 76

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