Recovering age-metallicity distributions from integrated spectra: validation with MUSE data of a nearby nuclear star cluster

Current instruments and spectral analysis programs are now able to decompose the integrated spectrum of a stellar system into distributions of ages and metallicities. The reliability of these methods have rarely been tested on nearby systems with resolved stellar ages and metallicities. Here we derive the age-metallicity distribution of M54, the nucleus of the Sagittarius dwarf spheroidal galaxy, from its integrated MUSE spectrum. We find a dominant old (8-14 Gyr), metal-poor (-1.5 dex) and a young (1 Gyr), metal-rich (+0.25 dex) component - consistent with the complex stellar populations measured from individual stars in the same MUSE data set. There is excellent agreement between the (mass-weighted) average age and metallicity of the resolved and integrated analyses. Differences are only 3% in age and 0.2 dex metallicitiy. By co-adding individual stars to create M54's integrated spectrum, we show that the recovered age-metallicity distribution is insensitive to the magnitude limit of the stars or the contribution of blue horizontal branch stars - even when including additional blue wavelength coverage from the WAGGS survey. However, we find that the brightest stars can induce the spurious recovery of an old ($>8$ Gyr), metal-rich (+0.25 dex) stellar population, which is otherwise not expected from our understanding of chemical enrichment in M54. The overall derived stellar mass-to-light ratio of M54 is M/L$_{\mathrm{V}}=1.46$ with a scatter of 0.22 across the field-of-view, which we attribute to the stochastic contribution of a young, metal-rich component. These findings provide strong evidence that complex stellar population distributions can be reliably recovered from integrated spectra of extragalactic systems.


INTRODUCTION
Analyzing ages and metallicities or other chemical abundances of any stellar systems give us insight into their assembly history. Galaxies and other stellar objects (e.g. nuclear star clusters (NSCs), see Neumayer et al. 2020, for a review) are assembled by a combination of in-situ secular processes such as star formation, as well as ex-situ accretion of other systems. The vary-ing relative contribution of these two processes will lead to complex stellar populations, thus detecting and quantifying their distribution in age and metallicity is crucial to understand how galaxies and other stellar systems assemble their stellar mass.
Deep color-magnitude diagrams (CMDs) still count as the most reliable and detailed view that we can obtain about stellar populations present in any stellar system. Spectroscopic follow-up studies of the individual stars can then also provide radial velocities and chemical abundances. Resolved CMD analysis however automatically restricts us to within the Local Group ( 1 Mpc), as greater distances make it impossible to resolve stars at arXiv:2005.03682v1 [astro-ph.GA] 7 May 2020 or below the main sequence turn-off, where most of the age information lies. As an example, the most detailed studies of stellar populations in NSCs are restricted to the nearby ones in the center of the Milky Way (MW: Do et al. 2009Do et al. , 2013Feldmeier et al. 2014;Feldmeier-Krause et al. 2015, 2017a and the Sagittarius dwarf spheroidal galaxy (M 54: Siegel et al. 2007;Bellazzini et al. 2008;Alfaro-Cuello et al. 2019).
Because integrated broad band colors are severely prone to age-metallicity-reddening degeneracies (e.g. Worthey 1994;Carter et al. 2009), integrated spectra present our most detailed view of unresolved stellar systems, as different absorption lines as well as the shape of the stellar continuum respond more sensitively to age and abundance pattern changes . With the advancement in sensitivity, wavelength coverage and spectral resolution of spectrographs (e.g. XShooter: Vernet et al. 2011) or the new generation of integral-field units (e.g. MUSE: Bacon et al. 2014), the analysis of integrated spectra is changing from line strength analysis (e.g. Worthey 1994;Worthey & Ottaviani 1997;Thomas et al. 2003Thomas et al. , 2005Trager et al. 2000a,b) to fitting the whole observed spectrum maximizing the information present in each spectral pixel (e.g. Cid Fernandes et al. 2013;Wilkinson et al. 2015;McDermid et al. 2015;Comparat et al. 2017;Goddard et al. 2017;Chauke et al. 2018;Kacharov et al. 2018).
Even though integrated line strength analysis is still crucial in understanding certain (galaxy) features sensitive to specific spectral lines such as the initial mass function (IMF) or individual α-element abundances (e.g. Martín-Navarro et al. 2018, it has been shown that single stellar population (SSP) equivalent ages and metallicities are biased towards the youngest stellar population present in the integrated light (e.g. Serra & Trager 2007;Trager & Somerville 2009). With full spectral fitting methods however, we are slowly moving towards uncovering the whole chemical enrichment history of a stellar system by fitting a linear combination of multiple SSP models to the observed spectrum. Yet, the techniques necessary to not only recover mean ages and metallicities, but also a distribution in that parameter space (e.g. Cappellari & Emsellem 2004;Cid Fernandes et al. 2005;Ocvirk et al. 2006a,b;Cappellari 2017;Wilkinson et al. 2017) are still under development. On top of that, additional factors influencing their performance, e.g. the wavelength coverage or spectral resolution are also not well understood yet. However, deriving age-metallicity distributions are crucial in deciphering the mass assembly of extragalactic systems (e.g. Boecker et al. 2019).
Given these challenges, it would be beneficial to study resolved and integrated spectra of the same system to test the reliability of age-metallicity distribution recovery with full spectral fitting methods. Here, we utilize the richness of the 3.5 × 3.5 MUSE (Multi-Unit Spectroscopic Explorer: Bacon et al. 2014) data set of M 54, the nucleus in the Sagittarius dwarf spheroidal galaxy, in order to apply the full spectral fitting to its integrated spectrum aiming to recover its multiple populations and compare the results to the resolved study of the same data set (Alfaro-Cuello et al. 2019). With respect to similar comparisons between resolved and integrated SFHs (Ruiz-Lara et al. 2015, our advantage lies in the power of the MUSE instrument (see also Kuncarayakti et al. 2016): a) the two approaches are performed using the same data, which means that possible instrumental effects stay the same, b) the metallicity of each star can be directly derived from its individual spectrum, which are regarded intrinsically more reliable than photometric metallicities and thus reduce degeneracies further, when stellar ages are determined from isochrone fitting. This paper is organized as follows: in Section 2 we present the three different integral field unit (IFU) data sets analyzed in this work and briefly describe the analysis of the resolved stars from Alfaro-Cuello et al. (2019); in Section 3 we describe our analysis method of deriving age-metallicity distributions from integrated spectra; in Section 4 we show the results of this technique in dependence of different integrated spectra of M 54; in Section 5 we compare our integrated analysis method with the resolved star analysis; in Section 6 we discuss our results and in Section 7 we give our conclusions of this comparison exercise.

IFU DATA OF M 54
In this Section we briefly describe the data we use to analyze the stellar populations of M 54 from its integrated spectrum. While our comparison between resolved and integrated stellar population extraction is focused on the MUSE WFM data from Alfaro-Cuello et al. (2019), we also include data from the MUSE WFM-AO science verification and the publicly available WAGGS (WiFeS Atlas of Galactic Globular cluster Spectra; Usher et al. 2017) survey exploiting different instrument systematics like wavelength coverage and spectral resolution. In Figure 1, we show the pointings on M 54 of the three data sets as an overview.  an extension of about 2.5 times the cluster's effective radius (r ef f = 0.82 = 6.78 pc (Harris 1996(Harris , 2010 at a distance of 28.4 kpc (Siegel et al. 2011)). MUSE is an integral field spectrograph, mounted on the UT4 of the Very Large Telescope at the Paranal Observatory in Chile. It has a wavelength coverage of 4750 − 9300Å with 1.25Å/pix sampling and a mean spectral resolution of ∼ 3000. More information about the observing strategy and data reduction can be found in Alfaro-Cuello et al. (2019).
To derive individual ages and metallicites for the member stars of M 54 and identify multiple populations, Alfaro-Cuello et al. (2019) performed five essential steps, which are stated here in short for clarity: 1. Extract individual spectra of the resolved stars with a wavelength dependent PSF-weighting technique (PampelMuse: Kamann et al. 2013) by using a photometric reference catalogue from HST (Siegel et al. 2007). All spectra with a SNR < 10 are excluded from further analysis. This corre-sponds to a limiting magnitude of I = 22 mag, which includes stars just below the turn-off region.
2. Fit all extracted spectra with a full spectral fitting software (ULySS: Koleva et al. 2009 5. Perform Gaussian mixture models in the derived age-metallicity parameter space to determine the most likely number of distinct stellar populations present in M 54 as well as the probability of each member star belonging to one of those populations.

MUSE WFM-AO
MUSE WFM is also offered with ground layer adaptive optics correction by the GALACSI (Ground Atmospheric Layer Adaptive Corrector for Spectroscopic Imaging) module aimed to double the ensquared energy in one 0.2 × 0.2 spaxel as compared to natural seeing. Due to the four laser guide stars the wavelength range around the NaD lines (5820 − 5970Å) is blocked. We include the analysis of MUSE WFM-AO science verification data (60.A-9181(A), PI: Alfaro-Cuello) of M 54 investigating, whether the Na notch filter has an influence on the stellar population inference. The data consists of a single, central 1 × 1 pointing.

WAGGS survey
MUSE does not cover blue wavelengths short of 4750 A, which is often raised as caveat (see planned Blue-MUSE: Richard et al. 2019), since it is commonly understood that young stellar populations (< 1 Gyr) or certain stellar evolutionary stages, such as horizontal branch stars, dominate at these bluer wavelengths (< 4000Å). It also misses other important spectral lines like the Ca H & K lines, at 3969Å and 3934Å respectively, which are particularly sensitive to metallicity changes at fixed temperature. Therefore their potential contribution to the overall integrated light might not be significant enough in the MUSE wavelength range. To test this, we additionally analyze the integrated spectrum of M 54 from the WAGGS survey (Usher et al. 2017).
The goal of the survey is to provide a library of globular cluster (GC) spectra in the Milky Way and its satellite galaxies with a higher resolution (R ∼ 7000) and wider wavelength coverage (3270 − 9050Å) than other studies to investigate their stellar populations in detail. For this purpose, they utilize WiFeS (Wide Field Spectrograph), a dual arm integral field spectrograph, on the Australian National University 2.3m telescope at the Siding Spring Observatory, which has a field of view of 38 × 25 , hence targeting the center of the GCs. The spectrograph offers four, high resolution gratings, U7000 (3270 − 4350Å with 0.27Å/pix), B7000 (4170 − 5540 A with 0.37Å/pix), R7000 (5280 − 7020Å with 0.44 A/pix) and I7000 (6800 − 9050Å with 0.57Å/pix) in order to achieve the large wavelength coverage. For more information about the observing strategy and data reduction see Usher et al. (2017). The integrated spectra for all their observed GCs are publicly available on their website 1 .
As the spectra from WAGGS consist of four parts corresponding to the four gratings, we determine the new flux in the overlapping regions as the error weighted mean in order to generate one continuous spectrum. The corresponding new inverse error in that region is the mean of the inverse of the two overlapping error spectra. Then we re-sample the entire spectrum to the highest pixel dispersion of 0.57Å/pix (using SpectRes: Carnall 2017).

METHOD FOR ANALYZING AN INTEGRATED SPECTRUM
In this Section, we briefly describe our approach of full spectral fitting and how multiple stellar populations in age-metallicity space can be derived from an integrated spectrum. In theory, there are only two ingredients needed: a single stellar population spectral library and a fitting machinery that fits the models to the data.

Full spectral fitting method
We can extract stellar populations properties from an integrated spectrum by viewing the integrated light as a linear combination of many single stellar populations, each with a different, single age and metallicity 2 . Hence, the full spectral fitting algorithm finds the optimal weight for each SSP spectrum, such that their sum best represents the observed integrated spectrum. Since SSP models are normally normalized to one solar mass, the best-fit weights are mass fractions.
Due to the vast parameter space of the SSPs (typically > 500 models) and the typical age-metallicity degeneracy, this inverse problem is usually ill-posed and ill-conditioned. In our case, ill-posed means that the solution is not necessarily unique, as many different SSP combinations can represent the data equally well. Illconditioned refers to the fact that fluctuations on the noise-level in the data can drastically change the solution. Regularized least-squares minimization is a common way to treat both of these issues. This technique is also implemented in the program pPXF (Cappellari & Emsellem 2004;Cappellari 2017), which we will use in this study. It has the advantage of being able to derive two dimensional age-metallicity distributions, which is not the case for other fitting algorithms that use similar approaches (e.g. STECKMAP: Ocvirk et al. 2006a,b).
In the context of pPXF, regularization provides a way to smoothly link the sparsely returned unregularized mass fractions in age-metallicity space until a certain criterion on the data fidelity is met. This smoothing may be well motivated in the case of galaxies, where we assume that chemical enrichment does not likely occur in a discrete manner. However, if the solution requires a more bursty star formation history, and the data quality is high enough, regularization will allow for that as well (Cappellari 2017).
How much the weights become smeared out is controlled by the regularization parameter (λ), whereas the way they are distributed is imposed by the regularization matrix (B), which is typically a finite difference operator. We note that we use the first finite difference (B = diag(1, −1)) throughout this work, however we made sure that the second and third order one give consistent results (see e.g. Ocvirk et al. 2006a;Huang et al. 2016;Kacharov et al. 2018;Boecker et al. 2019).
In any case, the essential pPXF fitting procedure is always the same, which follows the instruction in the source code of pPXF as well as Press (2007, Chapter 19.5). First, an unregularized fit is performed in order to re-scale the noise vector, such that this fit has a reduced χ 2 of unity. Then regularized fits are performed, tuning the regularization parameter such that the χ 2 of the regularized fit moves one standard deviation away, which corresponds to √ 2 · #pixel. This sets a regularization parameter, which allows for a maximum amount of smoothness that is still regarded to be consistent with the data (see e.g. McDermid et al. 2015;Kacharov et al. 2018;Boecker et al. 2019).
Before the first fit, we conduct the following, preparatory steps: SSP models are broadened to the wavelengthdependent line spread function (LSF) of MUSE 3 , as the SSP models have slightly higher spectral resolution in the blue part of the spectrum than MUSE. In case of the WAGGS spectrum, we broaden the observed spectrum to a constant FWHM of 2.5Å, which is the approximate average spectral resolution of MUSE. After that the wavelength range of interest is selected and residual skylines are masked out. Solely multiplicative polynomials are used to correct for any continuum mismatch between the SSP models and the observed spectrum. Their degree is determined according to (λ max − λ min )/200Å ensuring that spectral features narrower than 200Å are not influenced by the polynomial. For every fit we also include the corresponding error spectrum as determined by the data reduction process. In pPXF, all inputs have to be logarithmically re-binned in wavelength before the fitting process.
Together with the mass weights of the bestfit returned by pPXF and the predictions for the total luminosity of each SSP model in a certain photometric band, we can calculate the mass-to-light ratio (M/L) of the stellar system. We follow Cappellari et al. (2013, equation where, for the i-th SSP model, w are the weights returned by pPXF, M +rem is the mass in stars and dark remnants 4 and L V the corresponding V-band luminosity. We can also estimate the total stellar mass of the system from the integrated spectrum. For this we need to take into account the distance d to the stellar object, as well as normalization constants applied to the observed spectrum N obs and the SSP models N SSP prior to fitting 5 . Hence, we arrive at the following formula (see also Wilkinson et al. 2017;Kacharov et al. 2018) for the total stellar mass M , tot :

Single stellar population models
Due to the large age and metallicity as well as wavelength coverage (1680 − 50000Å) we chose the SSP models from the E-MILES library (Vazdekis et al. 2016) for our main analysis of M 54's stellar populations. Using the BaSTI isochrones (Pietrinferni et al. 2004) they cover 53 age bins between 0.03 − 14.0 Gyr and 12 metallicity [M/H] bins between -2.27 and 0.4 dex. We are using the bimodal IMF with a slope of 1. 3 Vazdekis et al. (19963 Vazdekis et al. ( , 2003, which is similar to a Kroupa (Kroupa 2001) IMF.
The spectral resolution of the E-MILES models is 2.51 A (FWHM) in the MUSE wavelength range until 8950.4 A, after that it jumps to about 4.2Å and then increases slowly with wavelength (see Vazdekis et al. 2016, Figure  8). This means that the LSF of MUSE is much narrower between 8950.4−9300Å than the SSP models and therefore we truncate the MUSE spectrum there. Additionally, the contamination from sky residuals can be quite significant in this regime anyway. How this choice of truncation influences the age-metallicity recovery is investigated in Appendix C.
In Section 6.1 we will discuss the impact of using different SSP models for our analysis. However, it is beyond the scope of this paper to provide a full comparison among different SSP libraries and their impact on the recovery of the age-metallicity distribution from integrated spectra. Thus, the reader is referred to more detailed studies of different spectral synthesis assumptions and techniques in the context of full spectral fitting (e.g. Conroy et al. 2009Fan et al. 2016;Baldwin et al. 2018;Dahmer-Hahn et al. 2018;Ge et al. 2018Ge et al. , 2019.

STELLAR POPULATION RESULTS FOR DIFFERENT INTEGRATED SPECTRA OF M 54
The flux-calibration of MUSE data makes it straightforward to construct an integrated spectrum of M 54. We can either sum up the individually extracted stars or collapse the whole data cube along the spatial axes.
In this Section we will probe the recovery of M 54's multiple stellar populations depending on the details of constructing the integrated spectrum from the MUSE data and on instrumental effects like wavelength coverage and spectral resolution by analyzing the WAGGS spectrum.
We summarize in Table 1 the different integrated spectra of M 54 analyzed in this work.

Integrated spectra from the entire MUSE cube
Firstly, we construct an integrated spectrum from the full cube for each of the 16 MUSE WFM pointings (see data set B in Table 1). Only spaxels with a formal SNR > 10 are included in the total integrated spectrum in order to avoid heavy sky residuals in the final integrated spectrum. Those are especially apparent for outer pointings. We do not exclude any possible contaminating sources and tested that the SNR cut does not affect the stellar populations recovery. We additionally combine these 16 integrated spectra into one single integrated spectrum (see data set A in Table 1). Note that we do not account for the overlapping regions of the 16 individual pointings (see Figure 1). The same is done for the single MUSE WFM-AO cube (see data set C in Table 1).
We then feed these spectra into pPXF with the E-MILES models considering all available ages and metallicities (636 models in total). The results from the pPXF fits are shown in Figure 2. We can see that the residuals are on the order of 2% emphasizing the excellent data and model quality.
Both integrated spectra from the 3.5 × 3.5 MUSE WFM mosaic and the single WFM-AO pointing are almost identical to the eye, and the recovered mass weights in age-metallicity space show a very similar distribution. A quantitative comparison between the recovered mass weights of both data sets is shown in Figure 11 a) in Appendix A.
We can identify an old, metal-poor (∼ 1.5 dex) stellar population at 8 and 14 Gyr, and a young (1 Gyr) and metal-rich (+0.25 dex) contributions. We also see a smaller contribution of old, but metal-rich mass weights. These weights do not fit into our astrophysical picture of chemical enrichment, hence their origin is further explored and discussed in Section 4.6.
Here, it is important to see that the lack of the sodium region in the WFM-AO data does not influence the recovery of the stellar populations properties. For the WFM data we also conducted tests of masking and not masking the NaD lines in the same region as Alfaro-Cuello et al. (2019), as this line is significantly broader than other lines in the spectrum due to interstellar absorption. Both showed consistent results meaning that pPXF is robust against such "outlier" spectral lines, at least if the wavelength range covers enough other prominent lines. The NaD lines can often be problematic as it is influenced by many different effects, such as the interstellar medium, IMF and sodium abundances.
We also show the recovered age-metallicity distribution from the analysis of the integrated spectrum for the 16 single pointings of the MUSE WFM data set in Figure 3. Overall, the three populations are picked up in every pointing with varying relative strength (see Figure 11 b) for a more quantitative comparison). The only significant outlier belongs to pointing number 12 (red color code), where the integrated spectrum is heavily influenced by a red supergiant star (V = 17.15 mag & I = 13.43 mag). It can also be clearly identified in Figure 1, as it is extremely red. This star contributes about 4% to the total flux of the entire MUSE pointing. If we mask this star out and re-do our integrated light analysis, we obtain a consistent age-metallicity distribution (blue color code) compared to the other pointings.

Integrated spectra from individual stars
The individually extracted stellar spectra from Alfaro-Cuello et al. (2019, see also Section 2.1) of the same 3.5 × 3.5 MUSE data set allow us to uniquely combine or exclude certain stars from the integrated spectrum and study the effect on the stellar population recovery. Only stellar spectra with SNR ≥ 10 are considered in a The SNR was estimated from the pPXF fit residuals between 5000Å and 5500Å. Due to the correlations in the noise and systematic effects from the data reduction, the signal-to-noise is heavily overestimated, when determined from the formal error cube. We therefore follow this more conservative approach as done in García-Benito et al. (2015); Sarzi et al. (2018). b The SNR for the U,R,B,I gratings respectively was taken from Usher et al. (2017, Table 2). Although the SNR seems to be higher than from the MUSE observations, by looking at the residuals of Figure 7 it becomes apparent that they are of the same order.
this sample and they make up roughly 50% of the total flux from the 3.5 × 3.5 MUSE mosaic.
First, we simply sum up the entire sample of extracted stars (7165 stars) and in a second test case we only consider those stars that were identified as member stars (6656 stars) (see data set C in Table 1). Identical stars that were extracted from different pointings were only accounted for once in the integrated spectrum by taking the corresponding stellar spectrum with the higher SNR. The flux of the individual stellar spectra is preserved and not normalized in any way prior to creating an integrated spectrum from them.
Results from the pPXF fit are shown in Figure 4. The recovered stellar populations in age-metallicity are nearly identical by eye, which is also shown in a quantitative comparison in Figure 11 a) in Appendix A. Hence, pPXF seems to be robust against contaminating sources, but the stars classified as non-members only make up 8% of the total integrated light. Still, these non-member stars can have radial velocities of up to -200 km/s, whereas the systemic velocity of M 54 is around 141 km/s (see Figure 3 of Alfaro-Cuello et al. 2019), which pPXF compensates with a larger velocity dispersion. The fitted velocity dispersion for the integrated spectrum including all extracted stars is around 15 km/s, whereas for the integrated spectrum only including member stars it is 1 km/s, which is the hard coded lower limit in pPXF. The true internal velocity dispersion of M54 is less than the spectral resolution of MUSE. Nevertheless, there is no apparent change in the recovered stellar population properties.

Integrated spectra with limiting magnitude cutoffs
We also look at different magnitude cuts in the CMD and its effect on the recovered stellar populations parameters, as it changes the relative contribution of stars in different evolutionary stages to the total integrated light. This is particularly interesting, if we think that certain regions of the CMD, like for example the horizontal branch, can be responsible for erroneously recovered stellar population properties. Even though making magnitude cuts is inconsistent with a SSP, we expect that the difference becomes negligible at a certain magnitude.
We show in Figure 5 three different cases, where we include all member stars with I-band magnitudes brighter than 16, 18 and 20 mag respectively (see data set E in Table 1). These cuts encompass 53%, 89% and 98% of the total integrated light of all member stars.
For the first case, a), only stars brighter than the red clump are considered, for which pPXF still finds the old and intermediate-age, metal-poor population, but no longer the young, metal-rich one. This is likely due to the fact, that the fraction of young, metal-rich stars contributing to the total spectrum is lower in this particular space of the CMD. As a consequence, there are more weights in the old (> 8 Gyr), metal-rich (> 0.0 dex) age-metallicity regime.
For the other two cases, b) & c), the result becomes identical to Figure 4, where all member stars were taken into account. This can also be seen in Figure 11 c) in Appendix A. In particular for case b), pPXF is able to reproduce the same results, as if the data had a limiting magnitude of 18 in the I-band, which is just below the red clump. This is not surprising, as the magnitude cut at 18 already includes 89% of the total integrated light coming from all member stars. due to interstellar absorption, however the age-metallicity recovery from pPXF is robust against the in-or exclusion of this region. The right panel shows the derived mass fractions in age-metallicity space that make up the bestfit from pPXF. b): The same as for the top panel but showing the integrated spectrum from the single pointing MUSE WFM-AO data. We see that the recovery is insensitive to the blocked region from the Na notch filter.
Moreover, it is reassuring to see that relatively larger contribution from the horizontal branch to the total integrated spectrum in b) does not artificially induce any additional young populations. This is investigated further in the following section.

A closer look at M 54's horizontal branch
M 54 has quite an extended horizontal branch (HB) with a ratio of 0.75 (Georgiev et al. 2009, a HB ratio of 1 and -1 means only blue or red HB stars respectively), which has been argued to bias age determinations by about 2 Gyr or more (e.g. Lee et al. 2000;Schiavon et al. 2004;Colucci et al. 2009;Georgiev et al. 2012) or to cause spurious young populations in integrated light analysis (e.g. Ocvirk 2010). In particular, metal-poor globular clusters have been found to exhibit bluer horizontal branches, which have strong Balmer lines and hence can mimic young main sequence stars. Generally, HB stars are difficult to model in SSP models, as there are higher order parameters determining their morphology (e.g. Lee et al. 1994;Gratton et al. 2010).
In principle, the extended HB in M 54 could be the reason that we find the metal-poor population (∼ -1.5 dex) at 8 Gyr even though we know that it is likely older than that. The advantage of our MUSE data set is that we can test this hypothesis by constructing an integrated spectrum without M 54's HB stars and fit it with pPXF. We exclude all stars that either have I < 18.5 mag and V − I < 0.7 mag or I > 18.5 mag and V − I < 0.3 mag. . Distribution of pPXF-recovered weights in age-metallicity space for the 16 pointings from the MUSE WFM dataset of M 54 (blue color code). They are arranged in the same order as they appear on the sky in Figure 1. The red color code (the scale is the same as the blue color code) shows the age-metallicity distribution recovered with the extremely red supergiant star contaminating the spectrum. Note that the apparent presence of an old and metal-rich component is likely related to some of the brightest AGB stars (see Figure 8 and Section 4.6).
As can be seen from Figure 6 the recovered stellar populations properties are identical to the fit, where all members stars were included in the integrated spectrum (see Figure 11 c). Hence, we can conclude that the HB is not responsible for potentially shifting the metal-poor population to 8 Gyr. In Section 5 we give another possible explanation for this being due to the large oxygen abundances in M 54.
On the other hand, if we fit the integrated spectrum made of only HB stars, we indeed recover very young (< 1 Gyr), metal-poor and metal-rich alike, populations and an old (∼ 10 Gyr), metal-poor population. We can even attribute the recovered old component to red (V − I > 0.43 mag) and the very young components to blue (V − I > 0.43 mag) HB stars.
Interestingly, pPXF always recovers a very small mass fraction (< 1%, which corresponds to < 10% in light) in the youngest age bin (0.03 Gyr) from fitting the E-MILES models to the MUSE integrated spectrum with and without HB stars alike. Hence, we cannot attribute this spurious weight in the youngest age bin to the pres-ence of blue HB stars in the integrated spectrum as was found in Ocvirk (2010) 6 .

Bluer wavelengths from WAGGS
The influence of hot stars, such a young stars or evolved horizontal branch stars, starts to become dominant in bluer wavelengths (< 4000Å) than the MUSE range. To be conclusive that M 54's extended HB does not influence the age-metallicity distribution recovery as found in the previous section, we here analyze M 54's integrated spectrum from the WAGGS survey until 3500 A 7 (see data set F in Table 1).
In Figure 7 a) we show the pPXF fits to the entire WAGGS spectrum and the integrated spectrum from the MUSE-AO observations. We chose the latter for the  comparison as it is centered on M 54 and hence could easily be cropped to the same field-of-view as WiFes in order to eliminate changes in the stellar population recovery induced by the differences in the spatial coverage. It is reassuring to see that both spectra acquired with completely different instruments are almost indistinguishable. The same applies to the residuals of pPXF, also when the WAGGS spectrum is fitted at native resolution (see Figure 15 in Appendix D). The recovered mass weights in age-metallicity space fitted to the whole wavelength range of WAGGS, only the wavelength covered by MUSE and to the actual MUSE-AO observations are plotted in Figure 7 b), c) and d) respectively (see also Figure 11 d) in Appendix A for a direct comparison). Including these bluer wavelengths does not recover any artificial young populations below 1 Gyr, therefore we conclude that the horizontal branch stars have no influence on the stellar population recovery from the integrated light in our analysis of M 54. Quite contrary, pPXF now puts all the mass weights at the 1 and 8 Gyr old population and the one at 14 Gyr vanishes. It reappears though, if the MUSE wavelength range is considered with the WAGGS spectrum. This could either imply that we loose the ability to recover very old ages, if we include these blue wavelengths or that the 14 Gyr population is not real/robust. However, it is more likely that this is associated with the overall difficulty to distinguish between SSP models with ages of 8 Gyr and above at fixed metallicity. Therefore, the "two" populations at 8 and 14 Gyr could also be just one. In Figure 15 in Appendix D we show M 54's integrated spectrum from WAGGS and MUSE fitted by the PEGASE-HR models, suggesting a more extended old, metal-poor population.

The influence of the brightest stars
The recovery of an old (8 − 14 Gyr) and metalrich (+0.25 dex) component in all of the above agemetallicity distributions does not necessarily fit into our astrophysical understanding of chemical enrichment in the presence of the three other populations. We therefore investigated, whether this component is real. From Figure 4 it seems that this component might arise from a few individual pointings. To explore this further, we repeated the exercise by only considering the classified member stars of M 54. We associate each member star of M 54 to the MUSE pointing it was extracted from, create an integrated spectrum and fit it with pPXF. For pointing number 12 we again excluded the red supergiant star as before.
Overall, the three populations, the old, metal-poor (∼ 1.5 dex) at 8 and 14 Gyr, and the young (1 Gyr), Grey shaded areas are masked out regions. The corresponding derived mass fractions in age-metallicity space that make up the bestfit from pPXF are also shown. Already for I ≤ 18 mag, the same distribution is found as if all member stars were included in the integrated spectrum.   metal-rich (+0.25 dex), are much more clearly recovered as seen in Figure 8 (see also Figure 11 e) for a more quantitative comparison). The old (8 − 14 Gyr), metalrich (+0.25 dex) component in pointing numbers 4, 5, 13 and 14 from Figure 3 vanishes. Looking at those fits again, it is evident that foreground stars with large offsets in their line-of-sight velocities are disturbing the integrated spectrum and the fit when summing up the full MUSE cubes. However, this is not the case for the central pointings (6, 7, 10 and 11) and number 16. In fact, only including M 54's member stars in the integrated spectra makes the old, metal-rich population much more prominent in the central four pointings. We speculate that this could be due to two effects, a population of highmetallicity MW foreground stars or perhaps thermally pulsing AGB stars. The former possibility was tested by running the Besançon Milky Way model (Robin et al. 2003)  To explore the second possibility, we looked at the color-magnitude diagram of the member stars per MUSE pointing (see Figure 12 in Appendix B). Especially in the four central pointings, brighter (I < 14 mag) and much redder (V − I 1.7 mag) stars are found, that could be thermally pulsing asymptotic giant branch (AGB) stars. These brightest stars contribute around 20% of flux to the total integrated spectrum and therefore have a non-negligible influence on the shape and spectral features of the integrated spectrum. Their red continuum shape and typical spectral features like prominent TiO bandheads can easily be mimicked by old, metal-rich stellar populations, if they are not properly accounted for in the SSP models (see e.g. Maraston 2005;Maraston et al. 2006, for the influence of AGB stars). If they are indeed the source of an old, metalrich component then it might also explain why the component was not recovered from the WAGGS data (see Section 4.5), as we expect that the influence of AGB stars becomes weaker at bluer wavelengths.
Excluding these stars from the integrated spectrum of the four central pointings and repeating our analysis, we obtain the age-metallicity distributions as shown in Figure 8. The contribution of the old, metal-rich component decreases significantly from 20−30% to 0.3−2%. Evidently, the old, metal-rich component is not recovered in a significant amount when considering the integrated spectrum made from the full MUSE cube in three out of the four central pointings (Figure 3). We attribute this to the fact that the full cubes contain enough flux from fainter, unresolved stars, such that the contribution of the brightest stars drop to around 10%.
For pointing number 16 a single star (marked by a grey circle in Figure 12) is responsible for the recovery of the old, metal-rich component. This star also appears in pointing number 15, but its presence in the integrated spectrum of this pointing does not cause the old, metalrich component to appear -probably because its flux contribution is only 4%, whereas in pointing number 16 it is 8%.
Counterintuitivly, the recovered age-metallicity distributions for pointings 4, 5 and 8 in Figure 8 do not show the old, metal-rich component even though their integrated spectra include stars with I > 14 mag that contribute around 20% to the total flux. Nevertheless, we confirm that excluding these brightest stars from the integrated spectra from Figures 4 b) and 5 b) & c), makes the old, metal-rich component completely vanish, while for Figures 5 a) and 6 the contribution significantly decreases. Instead, the relative contribution of the young, metal-rich component becomes stronger in all cases.

INTEGRATED VS. RESOLVED AGE-METALLICITY DISTRIBUTION RECOVERY
As we have seen in the previous section, the recovery of the mass distribution in age-metallicity space does not heavily depend on the exact approach used to construct the integrated spectrum of M 54 -however the contribution of the brightest stars is possibly responsible for the recovery of unphysical mass weights for this particular system (Section 4.6). Therefore, we restrict our comparison to the resolved stellar population analysis from Alfaro-Cuello et al. (2019) to the results from the integrated spectrum built from M 54's members excluding the brightest stars (see Section 4.6).
In Figure 9 a) we show the age-metallicity relation derived from the single stars binned to the same agemetallicity grid as the one we use in our integrated analysis to ease comparison. We show the stars that belong to M 54 in a red color code, where as stars that have been characterized as outliers in blue color-code, as quantified by Gaussian mixture models in Alfaro-Cuello et al. (2019). It is also important to note again that the horizontal branch stars are not included in this result.
In Figure 9 b) we show our result from the integrated light analysis. Now, the color code shows the absolute mass contained in each SSP bin instead of the relative mass fraction (see equation 2). We show the bestfit mass bins as well as the 84th percentile from randomly resampling the residuals from pPXF, adding them to the bestfit spectrum and re-fitting 100 times. We chose to keep the regularization parameter fixed to the value derived for the bestfit, as opposed to use no regularization at all. This allows for studying the effect on the smoothening of individual mass bins in age-metallicity space purely due to random variation in the fitted spectrum, and not the change of the regularization parameter. We also use the variation of the derived mass fraction from the 16 nearly independent MUSE pointings to quantify how much their absolute value as opposed to their smoothening across the age-metallicity plane. The mean relative differences between the mass fractions recovered for each pointing and the total of all member stars is around 28% and 38% for age and metallicity respectively (see Figure 11 e)). They are on the same order as the ages and metallicities measured from the individual stars from Alfaro-Cuello et al. (2019), which vary by 44% and 26% respectively across the pointings. They are likely induced by various factors such us differing SNR of the integrated spectrum or stochastic sampling of certain stars.
By focusing on the red color code of Figure 9, we see that the resolved star analysis shows a nicely rising chemical enrichment as a function of time. Most stars are roughly between 10 and 14 Gyr in age and -2.0 to -1.0 dex in metallicity. A more spread-out population of stars lie between 3 and 5 Gyr and around -0.25 dex in metallicity, whereas more stars seem to concentrate again at around 2 Gyr and solar metallicity.
The integrated light analysis on the other hand shows three more concentrated populations at 14 Gyr and -1.5 dex, around 8 Gyr and -1.25 dex and at 1 Gyr and super solar metallicity respectively. The returned chemical enrichment is also less continuously rising, but shows a rather flat enrichment from 14 to 8 Gyr and then jumps to +0.25 dex in metallicity at 1 Gyr. However, the separation between the 8 and 14 Gyr is likely not real, but either induced by the general poor age resolution at old ages, or by the quite complex individual element abundances of the Sagittarius nucleus. As shown by Carretta . The bins correspond to the age-metallicity grid adopted by Vazdekis et al. (2016) using the BaSTI isochrones (Pietrinferni et al. 2004). The red color code shows the stars that were characterized to M 54, whereas the blue color code shows outliers. The mean age and metallicity are 9.24 Gyr and -1.01 dex respectively. b): Age-metallicity distribution from fitting the full, integrated spectrum of M 54. Here, we show the result for the integrated light of individual member stars spectra excluding the brightest stars (see Section 4.6). The colorbar indicates now how much absolute stellar mass is contained in each SSP bin. The red color code shows the distribution of mass bins belonging to the bestfit solution, whereas the blue color code shows the extent of the 84th percentile mass bins derived from randomly re-sampling the residuals. The mass-weighted mean age and metallicity are 9.53 Gyr and -1.21 dex respectively.  Milone et al. 2011) and hence the solarscaled BaSTI isochrones are inconsistent with low metallicity MILES stars. Unfortunately, this likely does not explain why the metallicity derived for the metal-poor population in the integrated analysis is consistent with the resolved study and the young population is supersolar compared to solar values from the individual M 54 stars. If α-abundances are causing the metallicities differences between the two methods, we would expect to observe the opposite trend. Furthermore, the fit residuals of the Mgb lines (see e.g. Figure 7) show that the magnesium abundance is actually overpredicted in the SSP models suggesting discrepancies between the adopted and actual alpha-abundance of M 54.
Furthermore, the discontinuity between the old, metal-poor and the young, metal-rich population in the integrated analysis can have a number of possible explanations. It could be that the relative contribution of stars between 2 and 8 Gyr as seen in the resolved analysis is not significant in the integrated spectrum to be picked up by pPXF. This means that the star formation rate of a certain star forming episode has to reach a specific threshold to be contributing significantly to the integrated light. It could also be that this is an issue of how regularization is applied, as it smooths the mass weights only in the horizontal and vertical direction in the age-metallicity plane, but not diagonally.
Despite the apparent differences (or similarities) between the age-metallicity distributions of the resolved and integrated light analysis, the average quantities of both methods agree well. We quote a mean age and metallicity of 9.24 Gyr and -1.01 dex for the resolved stars and a mass-weighted mean age and metallicity of 9.53 Gyr and -1.21 dex for the integrated method. This corresponds to a difference of only 3% in age and 0.2 dex in metallicity, which is good precision considering the range of metallicities of almost 2 orders of magnitude. Weighting the resolved stars by their V-band luminosity yields a mean age of 9.69 Gyr and metallicity of -1.11 dex, whereas light-weighted quantities from the integrated method produce a mean age of 7.34 Gyr and metallicity of -0.99 dex.
The averages across the 16 MUSE pointings vary by about 0.43 Gyr and 0.06 dex for the resolved and by 0.82 Gyr and 0.10 dex for the integrated analysis. Statistical errors on these quantities from both Alfaro-Cuello et al.
(2019) and our random re-sampling of the residuals are below a few percent. However, based on the general uncertainty of stellar population synthesis as well as the poor age resolution at old ( 8 Gyr) ages, we do not claim to recover the true mean age and metallicity of M 54 to better than 20%.
To summarize, the integrated light analysis of M 54 can recover a young, metal-rich and old, metal-poor stellar population even though pPXF is free to choose any, not necessarily physical, age-metallicity combination that best represents the observed integrated spectrum. The derived mean ages and metallicities are consistent with the resolved analysis despite the differences in the used SSP models and the lack of considering individual element abundances that are present in M 54.

Mass-to-light ratios from integrated analysis
From the returned pPXF mass weights we also calculate stellar M/L ratios in the V-band (see equation 1) for integrated spectra made from M 54 member stars (excluding the brightest stars; see Section 4.6). This is done for our canonical E-MILES SSP library choice with a bimodal IMF of slope 1.3. We find a global value M/L V = 1.46, however across the 16 MUSE fields we find values that vary from 1.3 to 1.8 (or with a standard deviation of 0.22 from the global M/L). These are shown in Figure 10 as a function of luminosity density colorcoded by contribution of the young (1 Gyr), metal-rich (+0.25 dex) component. From this it appears that the four central pointings, which have the highest luminosity density, tend to have lower M/L ratios. The outer fields at lower luminosity density show a wide spread in M/L ratios. We confirm that the scatter in the derived M/L ratios originate from the varying relative contribution of the young, metal-rich mass fractions, as this directly translated to a change in the M/L ratio (see equation 1). This agrees with the young, metal-rich stars being more centrally concentrated as found by Alfaro-Cuello et al. (2019). There is an outlier corresponding to pointing number 4, which has a low M/L (≈ 1.3) at low luminosity density (≈ 20 L /pc 2 ), likely caused by the brightest star contributing around 25% to the total flux (see Figure 12).
The exclusion of the horizontal branch in the composite spectrum (see Section 4.4 after 4.6) yields a M/L ratio of 1.45, whereas the magnitude cuts (see Section 4.3 after 4.6) yield 1.53, 1.38 and 1.48 for I < 16, 18 and 20 mag respectively. Even though these changes are below our estimated statistical errors of 1 − 2%, we do not claim that they are significant, especially with respect to typical uncertainties of 6% in other studies (e.g. Cappellari et al. 2013).
Our stellar population derived M/L V ratio agrees with measurements from Kimmig et al. (2015), who modelled M 54's velocity dispersion profile with a King profile and taking internal rotation into account. However, it is lower by about 0.5 compared to studies from Baumgardt & Hilker (2018); Dalgleish et al. (2020), who fitted N-body simulations without internal rotation to the velocity dispersion profiles. This might be an indication that it is important to include internal rotation in the dynamical modelling, as it decreases the M/L ratio measurement. There is evidence after all that M 54's young, metal-rich population exhibits a significant rotation signature (Bellazzini et al. 2008;Alfaro-Cuello et al. 2020).

DISCUSSING COMPARISONS OF INTEGRATED AND RESOLVED STUDIES
After thoroughly having compared the stellar population results from resolved star analysis of Alfaro-Cuello et al. (2019) and our integrated light analysis in the previous section, we have now arrived at the question, whether this comparison proves that the recovery of agemetallicity distributions from integrated spectra provides us with the same information content as the resolved stars with regard to stellar population ages and metallicities.
Assuming now that the derived ages and metallicities of the resolved stars resemble the 'best' knowledge we have about M 54, one might take the mismatch between the two approaches in Figure 9 as a failure of the integrated method. However, the importance and success of this comparison is not to be measured in how perfectly the individual bins in age and metallicity match each other, but the fact that the integrated light analysis can clearly and robustly detect that M 54 hosts multiple populations, which are even located in a similar agemetallicity space as the properties of the resolved stars. Considering that the resolved and integrated analysis  Figure 10. Mass-to-light ratios in the V-band derived from returned pPXF mass weights for integrated spectra made from M 54's member stars excluding the brightest stars (see Section 4.6) plotted against the luminosity density in each corresponding pointing. The symbol size corresponds to the distance from M 54's center, while the color-code shows the fraction of the young (1 Gyr), metal-rich (+0.25 dex) component. The orange star corresponds to the global M/L ratio obtained from an integrated spectrum made from all the M 54's member stars excluding the brightest stars. Estimated errors from randomly re-sampling the residuals is on the order of the symbol size.
techniques are also very different conceptually and in the models they use, the similarities of the recovered age-metallicity distributions are compelling.
In the following sections we will discuss in more detail, which aspects of stellar population analysis have to be considered to perform a one-to-one comparison between the resolved and integrated methods (Section 6.1). Furthermore, we argue that neither of those two techniques should be regarded as more reliable than the other (Section 6.2).

Is the comparison actually self-consistent?
Our comparison of integrated light versus resolved stellar population studies has the main advantage that it uses the same dataset. Nevertheless, we still need to consider two aspects related to the fitting of the stellar populations, that are both nontrivial to implement: First, we would need to make sure to use the exact same stellar synthesis models. Only then, we would be able to estimate, if discrepancies in the derived age and metallicity properties between the two techniques are induced by the different models or the different approaches themselves. The resolved study from Alfaro-Cuello et al. (2019) uses the ELODIE 3.2 (Prugniel & Soubiran 2001Wu et al. 2011) stellar library to estimate the stellar parameters for 4750 − 6800Å and the Dartmouth (Dotter et al. 2008) isochrones in order to determine the stellar ages. On the other hand, we have used the E-MILES SSP library (Vazdekis et al. 2016) together with the BaSTI isochrones (Pietrinferni et al. 2004) to fit for ages and metallicities in the wavelength range of 4750 − 8950.4Å. In this wavelength regime, the E-MILES SSP models are based on three different stellar libraries: the MILES (Sánchez-Blázquez et al. 2006), the near-IR CaT (Cenarro et al. 2001) and the Indo-U.S. (Valdes et al. 2004) library (see also Vazdekis et al. 2012). With these differences basically all systematics regarding spectral and stellar synthesis as well as their respective modelling are captured. Different sets of isochrones have different assumptions about stellar evolution, in-or exclude certain evolutionary phases of stars and are computed for different age and metallicity bins. Different stellar libraries have different flux calibration, wavelength coverage and, in case of empirical ones, are biased towards metallicities in the solar neighbourhood. All of this can influence the derived absolute ages and metallicities of stellar populations in both analysis techniques neglecting additional systematics induced by, for example, individual element abundances (Vazdekis et al. 2001;Schiavon et al. 2002) as discussed in Section 5.
In principle, this issue could be resolved by using the exact same models in both approaches. Where the authors of the resolved analysis have in theory the freedom to choose any combination of stellar libraries and set of isochrones, analysis of integrated light is restricted to the publicly available SSP models, which have a fixed combination of stellar library and isochrones 9 . Changing models and re-doing the integrated analysis is straightforward, but in the resolved case steps 2. − 5. in Section 2.1 have to be repeated, which can be quite time consuming for a decent amount of models and are out of the scope of this study.
For the interested reader, we show in Appendix D our integrated analysis of M 54 conducted with the PEGASE-HR SSP library (Le Borgne et al. 2004), in order to match the stellar library used in resolved method, however there the isochrones are from PADOVA (Bertelli et al. 1994). It is quite interesting to see how much the recovered age-metallicity distribution seems to depend on the adopted SSP models at first glance, while they are still recovering the same physical implications of M 54's multiple populations, all modelling uncertainties considered.
The second aspect that arises when trying to compare the two approaches is that the color code in Figure  9 does not represent the same physical quantity. The resolved study shows the number of stars in each agemetallicity bin, but the integrated analysis provides us with a mass (or light) fraction of an SSP corresponding to a particular age and metallicity. From the latter, we can deduce the absolute mass in each bin and consequently the total mass of M 54 relatively easily. Inferring the mass of each star can be deduced from the fitted isochrone and their magnitudes (see e.g. also Pont & Eyer 2004;Lin et al. 2018). However, a completeness correction is necessary to account for non-detected stars below the turn-off, where most of the mass lies. Only then could the stellar mass in each age-metallicity bin for the individual star analysis be estimated. The individual star counts also influence the derived mean ages and metallicities making the comparison of their values to the integrated measurements not one-to-one let alone in a light-or mass-weighted sense.
6.2. Is one method more reliable than the other?
After discussing the potential ways of making the comparison between the resolved and integrated analysis of M 54 as self-consistent as possible, the question arises, whether we would gain new knowledge from this extra work. Certainly, this would put the two approaches to the ultimate test, but we would also need to assume that one approach is better or more reliable than the other. In fact, both approaches suffer from the same difficulties that are connected to the well-known degeneracies and difficulties in stellar population modelling, such as the age-metallicity degeneracy, unknown individual element abundances and poor age resolution at old ages ( 8 Gyr).
With regard to the resolved stars from the MUSE data, their iron abundance is generally well-determined by their spectra, however the age determination from isochrone fitting is rather degenerate, as especially other element abundances shift the isochrones and can induce artificial age variations (see Section 5). Furthermore, fitting an isochrone through one point on the CMD is very degenerate in itself, as can be seen from the outliers in Figure 9 a). Some kind of measure needs to be defined in order to identify these outliers, as was done in Alfaro-Cuello et al. (2019) with the means of Gaussian mixture models.
Similarly, in the integrated analysis, SSP models at old ages ( 8 Gyr) and at fixed metallicity are more or less indistinguishable, therefore the age leverage is poor in this regime. This essentially could mean that the two old metal-poor populations in Figure 9 b) (one at around 8 Gyr and the other one at 14 Gyr) are the same, which is further complicated by the high oxygen abundances in M 54 as discussed in Section 5. Furthermore, here it is also hard to tell, which mass weights are robust and which are erroneously being generated simply due to the ill-posed nature of the inversion problem. Nevertheless, the metallicity determination seems to be more robust and hence having a handle on the metallicity distributions from integrated spectra of extragalactic objects is already a big advantage as compared to average values.
A similar argument holds when it comes to the comparison of different SSP models with the same technique. We do not know which models are intrinsically more reliable or closer to truth in nature, although a good parameter coverage across log g, T eff and [Fe/H] (and potentially [α/Fe]) of stellar spectra is always a limiting factor. A choice of a certain SSP library always has to be made, therefore absolute quantities, such as the exact position of the mass weights in age-metallicity space might not be as reliable. However, the relative trends are expected to be the same. Meaning that, if we always use the same SSP models for different stellar objects of interest, we will be able to say differentially, whether the objects have experienced different chemical enrichment histories.
In conclusion, without being prejudiced against the credibility of either of the two methods, we can reliably say that both results in Figure 9 show the following results: 1. There are multiple stellar populations present in M 54.
3. The cluster is dominated by the old, metal-poor population.
4. The mean age and metallicity are in the range of 9 − 9.5 Gyr and -1.0 to -1.2 dex.
Hence, the integrated analysis is capable of identifying multiple stellar populations from a single integrated spectrum. It results in a similar star formation and enrichment history as the resolved analysis based on CMD analysis and spectral fitting of individual stars.

CONCLUSIONS
In this work we have presented the analysis of M 54's integrated spectrum from three different data sets (MUSE WFM, MUSE WFM-AO and WiFes) with the goal to recover its multiple stellar population content via full spectral fitting (pPXF: Cappellari & Emsellem 2004;Cappellari 2017) of a library of SSP models (E-MILES: Vazdekis et al. 2016) to the observed spectrum.
Thanks to the individually extracted stellar spectra of the 3.5 × 3.5 MUSE WFM data set from Alfaro-Cuello et al. (2019), we could also investigate the influence on the stellar population recovery by excluding the contribution of certain stars to the total integrated spectrum. In light of all our tests, we draw the following conclusions in recovering age-metallicity distributions from integrated spectra: • The derived mass fractions in age-metallicity space are robust against 1) the Na notch filter in MUSE-AO observations (Figure 2), 2) the inclusion of stars classified as non-members ( Figure 4) and 3) the contribution of extended horizontal branch stars ( Figure 6).
• The recovery of the age-metallicity distribution is not very sensitive to the limiting magnitude of the observations. Consistent results are achieved, even if the limiting magnitude were 4 times brighter than the main sequence turn-off region ( Figure 5).
• The recovered mass fractions are consistent in their absolute position in age-metallicity space over individual pointings of the 4 × 4 MUSE mosaic, as long as the spectrum of an overly bright star does not dominate the integrated spectrum ( Figure 3).
• Bright (I < 14 mag) and red (V − I > 1.7 mag) stars in the integrated spectrum seem to induce erroneous old (> 8 Gyr), metal-rich (+0.25 dex) populations in the recovered age-metallicity distribution ( Figure 4). Uncertain evolutionary phases such as the thermally pulsing AGB not included in the SSP models could be an explanation for this.
• The absolute derived ages and metallicities change, as expected, with different SSP model assumptions, however differentially the trends stay the same (Figure 9, 14). Hence, M 54's multiple stellar populations are indeed retrievable from its integrated spectrum showing an old (8 − 14 Gyr), metal-poor (-1.5 dex) as well as a young (1 Gyr) and metal-rich (+0.25 dex) population.
• The derived mass-weighted mean age and metallicity of 9.53 Gyr and -1.21 dex are consistent with the corresponding averages of the resolved analysis of 9.24 Gyr and -1.01 dex respectively.
• The derived stellar M/L ratios show more stochasticity in the outer regions of M 54 (M/L V = 1.3 − 1.8), where the luminosity density is lower, as compared to the central region, where the value converges to around 1.46. We attribute this to the lower relative contribution of young, metal-rich mass fractions.
In this context we also compared and discussed our results with findings of the resolved stellar population analysis from same MUSE WFM data set (Alfaro-Cuello et al. 2019). From this we find that age-metallicity distributions can be derived from full spectral fitting of integrated spectra with comparable reliability as from resolved studies, as both approaches suffer equally from the same difficulties, uncertainties and degeneracies in stellar population synthesis modelling, especially with regard to age determinations, whereas the recovered metallicity distribution seems to be more robust. While IFU observations of resolved systems can certainly provide detailed information on a star-by-star basis, our integrated approach can provide the same information content, if the scientific goal is to disentangle multiple or complex stellar populations of stellar systems. It is also worth noting that the integrated analysis reveals this information with a single fit in several minutes as opposed to lengthy data extraction and analysis steps undertaken in the case of the resolved study (see Section 2.1). On top of that, our returned distributions have a physical unit attached to them (mass fractions) instead of number counts, which lets us straightforwardly calculate the stellar mass of the different populations or the system in total as well as M/L ratios. This provides us with a quick and detailed knowledge about the stellar content of an object. In spite of the modelling differences between both methods, we find that the average age and metallicity from the integrated and resolved stars analysis agree remarkably well with each other. This is of key importance for extragalactic studies at low and high redshift, which can only access the integrated light -especially now with the advanced development of chemo-dynamical models for external galaxies (e.g. Poci et al. 2019;Zhu et al. 2020).
With the ability to study the two dimensional distribution in the age-metallicity plane of thousands of stellar objects, we can establish a connection between the properties of multiple stellar populations to their global properties like total mass, presence of super-massive black holes and environment. This means, we are now in an era (data-and modelling-wise) to constrain formation scenarios of nuclear star clusters or the stellar mass assembly of galaxies on a statistical scale with the stellar population distributions from integrated spectra.
We would like to thank the anonymous referee, who provided a very careful and constructive revision that improved this manuscript. AB likes to thank Iskren Georgiev and Nikolay Kacharov for their help in this project. We also thank Anil Seth for useful discussions. This work was funded by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) -Project-ID 138713538 -SFB 881 ("The Milky Way System", subproject B08). IMN acknowledges support from the AYA2016-77237-C3-1-P grant from the Spanish Ministry of Economy and Competitiveness (MINECO) and from the Marie Sk lodowska-Curie Individual SPanD Fellowship 702607. RL also acknowledges support from the Natural Sciences and Engineering Research Council of Canada PDF award and DAAD PPP exchange program Projekt-ID 57316058. This research has made use of the services of the ESO Science Archive Facility.

A. QUANTITATIVE COMPARISON BETWEEN RECOVERED MASS FRACTIONS
We show in Figure 11 the results from fitting the different integrated spectra of M 54 from Section 4 in the form of one dimensional distributions as a function of age and metallicity respectively. This allows for a more quantitative comparison between the different fits than the two dimensional age-metallicity distributions. Evidently the recovered mass fractions show overall consistent results between the various integrated spectra that we investigated. Comparing panels b) and e) in Figure 11 we see that the derived mass fractions across all 16 independent MUSE pointings show less variations, especially in the young (1 Gyr), metal-rich (+0.25 dex) component, when foreground stars and the brightest member stars are not included in the resulting integrated spectrum.  Figure 12 shows the CMDs of M 54 member stars extracted from their corresponding MUSE pointings from the 3.5 × 3.5 mosaic (Alfaro-Cuello et al. 2019). Especially, the central pointings (number 6, 7, 10 and 11) have a high number of bright (I < 14 mag) and extremely red (V − I 1.7 mag) stars that contribute around 20% to the total flux of the integrated spectrum made from member stars in that corresponding MUSE pointing. When these stars are excluded from the integrated spectrum, the contribution of the previously recovered old (> 8 Gyr) and metal-rich (+0.25 dex) mass fractions decreases from 20 − 30% to under 2%. In pointing number 16 a single star with 8% flux contribution was responsible for pPXF to recover these mass weights. The same star is also present in pointing number 15, where it did not cause the old, metal-rich component to be picked up, possibly because its flux contribution is decreased to only 4%. Interestingly, other pointings (e.g. 4, 5 and 8) also have a few bright stars that contribute a significant amount to the total flux, but did not cause any old, metal-rich component to appear in the derived age-metallicity distribution. A possible explanation might be that these stars (expect the bright star in pointing 4) do not have red continuum shapes as well as TiO absorption bands that become visible in the integrated spectrum, which could easily mimic old, metal-rich stellar populations in the SSP models. As mentioned in Section 3.2 we here present the mass distribution in age-metallicity space derived from the integrated spectrum of the entire 3.5 × 3.5 MUSE mosaic and made from the individually extracted member stars (Alfaro-Cuello et al. 2019) fitted to the entire wavelength range of MUSE with the E-MILES models. The results are shown in Figure  13. We see that pPXF now recovers a third component, which is at 14 Gyr and the lowest metallicity of -2.27 dex. It is very prominent, when the integrated spectrum was made from all 16 MUSE data cubes combined. Apparently, the recovery of this component depends on the wavelength range between 8950.4 − 9300Å, which is quite noisy, as can be seen from the residuals. Moreover, in this wavelength regime the E-MILES models have a lower resolution (FWHM ≈ 4.2Å) than the MUSE spectrum, which makes the two Paschen lines (n = 9 and n = 10) appear very broad. Nevertheless, in the observed spectrum they do not appear nearly as deep as the bestfit model.
Following the discussion in Section 5 we argue that this does not have any influence on our statements regarding the ability and reliability of recovering multiple populations from integrated spectra, as the exact absolute values of the ages and metallicites depend on the adopted SSP models. We can still make the same qualitative statement about M 54's multiple stellar populations and now the overall trend of the recovered chemical enrichment is even more consistent with a steady rise than when the wavelength rage was cut off.
To assess whether the broader spectral resolution of the E-MILES in that particular wavelength range could cause the new stellar population component, we have convolved the integrated spectrum to the lowest spectral resolution present in the E-MILES library at those wavelengths (FWHM 4.4Å). Still the 14 Gyr old component in the lowest metallicity bin was recovered. Here, we provide results of fitting the integrated spectrum of M 54 with the PEGASE-HR models (Le Borgne et al. 2004). They are based on the ELODIE 3.1 (Prugniel & Soubiran 2001Prugniel et al. 2007) high resolution spectra (R = 10000) and the Padova isochrones (Bertelli et al. 1994). The wavelength coverage is 3900 − 6800Å, ages and metallicity span 0.001 − 20 Gyr (68 bins) and -2.3 to 0.69 dex (7 bins) respectively. The assumed IMF is Kroupa (Kroupa 2001) and the mass of one SSP is also normalized to unity. We set the minimum age bin to 0.1 Gyr and the maximum to 14 Gyr in order to match the boundaries of the E-MILES models, however we made sure that no spurious mass weights were detected when including the full age range in PEGASE-HR models. Similarly to Kacharov et al. (2018), we also performed a pPXF fit with PEGASE-HR models that were interpolated to a finer age-metallicity grid. The models were fit to the integrated spectrum made from the individually extracted member star spectra. Both results and the comparison to the resolved study from Alfaro-Cuello et al. (2019) are shown in Figure 14.
We again detect an old, metal-poor population with a metallicity between -1.0 and -1.5 dex and an essentially unconstrained age between 8 and 14 Gyr. Now, we can also identify an intermediate population at around 3 Gyr and between -1.0 and -0.5 dex in metallicity, which is still more metal-poor than the intermediate-age population from Alfaro-Cuello et al. (2019). We also still see the young, metal-rich population at 1 Gyr and around +0.25 dex, again more metal-rich than the one identified in the resolved analysis. Mass weights with this same super solar metallicity but older ages are again attributed to the same systematics as discussed in Section 5.
The result from the interpolated PEGASE-HR models shows mass weights that are in the same location in the age-metallicity space as the fiducial models, but are on average lower. The weights are hence smeared out across several more age-metallicity bins, which is not surprising as the interpolation adds more linearly dependent models into the design matrix.
In Figure 15 we show the results of fitting the WAGGS spectrum of M 54 with the fiducial PEGSASE-HR SSP models, once for the native WAGGS spectral resolution (FWHM 0.8Å) and once broadened to 2.5Å to mimic the spectral resolution of MUSE. Both results are almost identical, whereas the high resolution fit retrieves much less mass weights at old age and super solar metallicity and gives more weight to the old, metal-rich (around 10 Gyr and -1.5 dex) and the young, metal-rich (around 1 Gyr and +0.25 dex) population as compared to the result from the lower resolution spectrum. We can also see some differences in the distribution of the recovered mass fractions in age-metallicity space, if we compare these to the results for the MUSE spectrum fitted with the PEGASE-HR models in Figure 14. Nevertheless, this gives us confidence that the recovery of ages and metallicities of multiple stellar populations from integrated spectra is not severely dependent on the spectral resolution.
Even though the fit to the WAGGS spectrum with the PEGASE-HR models in Figure 15 includes bluer wavelengths than the MUSE data, we still recover an extended old, metal-poor population, where compared to fitting the full WAGGS wavelength range with the E-MILES models in Figure 7, the metal-poor mass fractions at 14 Gyr vanished and instead concentrated around 8 Gyr. We argue that this difference arises because of the diverse modelling assumptions of stellar population synthesis and not because of the inclusion of bluer wavelengths. The integrated analysis has been conducted with the fiducial PEGASE-HR models. c): The same comparison as b), but here the PegaseHR models where interpolated onto a finer agemetallicity grid prior to fitting. Here, the individual mass weights returned by pPXF are on average much lower than with the fiducial model.  The fit to the native spectral resolution (FWHM 0.8Å) is shown in blue, whereas the broadened one (FWHM 2.5Å) is in orange. The residuals are on the 2% level for both resolutions. b): Recovered age-metallicity distribution from the high resolution spectrum. c): Same as b), but for the broadened spectrum with a FWHM comparable to MUSE.