Far-ultraviolet to FIR Spectral-energy Distribution Modeling of the Stellar Formation History of the M31 Bulge (2024)

Table 1.Overview of Data Used

Download table as: ASCIITypeset image

The UVIT instrument onboard AstroSat has spatial resolution of ≃1'', a 28 arcminute field of view, and is capable of observing in a variety of FUV and NUV filters (Tandon et al. 2017). New in-orbit calibrations of UVIT were carried out by Tandon et al. (2020) and are applied here. The data processing was carried out using CCDLab (Postma & Leahy 2017, 2021), with updated astrometry calibration from Postma & Leahy (2020), and updated photometry from Leahy et al. (2021a, 2021b).

The M31 survey covers the sky area of M31 with pointings labeled Field 1 through 19 (Leahy et al. 2020). The M31 bulge is located in the center of Field 1. It was observed by UVIT in 2017 (labeled obervation A) in the N279N, N219M, F172M, and F148W filters. ¹ A three-color FUV–NUV image of Field 1 is presented as Figure 2 of Leahy et al. (2018). Field 1 was observed again in 2019 in the UVIT filters F172M, F169M, and F148W. The 2017 observations (labeled A) included N279N, N219M, F172M, and F148W filters, and the 2019 observations (labeled B) included F172M, F169M, and F148W filters. Merged images were created by adding the counts images for A and B, adding the exposure maps for A and B, then dividing the counts images by the exposure images to create a count/s image. The photometry here was carried out using the merged UVIT images. ²

The merged F148W filter image of Field 1, which has the longest exposure and highest signal to noise, is shown in Figure 1. The bulge of M31 is approximately the central 7 arcminute diameter part of Field 1. Beyond this central part, light from the spiral arms of M31 begins to dominate over the light from the bulge. Thus in order to minimize inclusion of light from the spiral arms we use a maximum radius (major axis) from the center of the bulge of 450 UVIT pixels, which is 18756.

From Leahy et al. (2021c), the bulge at FUV and NUV wavelengths is elliptical with minor to major axis ratio of ≃0.75 at position angle of ≃42° east from north. We use an ellipse with this axis ratio and position angle, shown by the outer ellipse marked in Figure 2. This ellipse was subdivided in order to study spatial variations by dividing into 10 concentric elliptical annuli, each with equal area (so the radius intervals are not equally spaced), then further dividing these annuli into four quadrants, with the division between quadrants oriented along the major and minor axes of the bulge, as shown in Figure 2.

In addition to the data for the five UVIT filters, we obtained archival image data for M31 in five filters from the Sloan Digital Sky Survey (SDSS), two filters from the Spitzer Space Telescope IRAC camera, three filters from the Herschel Space Observatory's SPIRE instrument, and two from the Herschel Space Observatory's PACS instrument. We used SDSS Data Release 16 (Ahumada et al. 2020). The SDSS-IV Overview is given in Blanton et al. (2017). The Sloan Foundation 2.5 m Telescope description is given by Gunn et al. (2006) and the SDSS photometry is described in Doi et al. (2010). The Spitzer Space Telescope is summarized in Werner et al. (2004) and the IRAC camera is described in Fazio et al. (2004). The Herschel data was obtained from the NASA/IPAC Infrared Science Archive. The Herschel Space Observatory is described in Pilbratt et al. (2010), and the SPIRE and PACS instruments are described in Griffin et al. (2010) and Poglitsch et al. (2010), respectively.

For each of the 17 different filter images listed in Table 1, data were extracted from 10 equal area elliptical annuli around the nucleus. The annuli had an axis ratio of b/a = 0.75 and a major axis position angle of 42° east of north. Because the edges of the SDSS image fields were offset just NE of the nucleus (bulge center), the common center for the annuli was shifted slightly NE so that each quarter was not too close to the edges of the SDSS image fields. ³

The sum of image values in each annulus in each image was extracted using SAOImage DS9, then converted to millijanskys (mJy) using the respective flux conversions for the appropriate filter bands for UVIT, SDSS, IRAC, and SPIRE. For PACS, the images were in units of Jy/pixel so no further conversion was necessary. Errors in the flux values are calculated as the square root of the instrument counts, converted back into flux. Additional errors are added in quadrature. The additional errors are as follows. For SDSS, the u filter has an additional 0.1 DN/pixel error due to a bias in the sky-level determination. For IRAC, the errors are ≃10% of image flux. For UVIT, the errors in the magnitude-flux conversion are listed in the instrument calibration paper (Tandon et al. 2020). For SPIRE/PACS, the image files include error images in the same units, thus errors per pixel are extracted for the annuli in the same way as the image, and converted to errors in the sum by dividing by square root of number of pixels in the region.

A new SED fitting program was written that takes as input a dense grid of age and metallicity values, with range of ages 3.98 × 10⁶ yr to 1.5 × 10¹⁰ yr and range of metallicities Z = 5.0 × 10⁻⁴ to 0.04162 (log(Z/Z_⊙) = −2.617 to 0.303, when taking Z_⊙ = 0.0207). The program carries out a two-dimensional interpolation on the grid to calculate the fluxes for a SSP with any specified age and log(Z/Z_⊙). A χ² minimization can be carried out for the cases of a single SSP (with 4 parameters: age, log(Z/Z_⊙), mass and extinction E(B-V)), two SSPs (with eight parameters, four for each SSP), or three SSPs. The program also carries out a parameter grid search around the best-fit set of parameters to estimate the uncertainties in the parameters.

An SSP model consists of a set of stars of a single age and metallicity. More realistic models are constructed by combinations of SSPs with different ages and metallicities. Here, we used the SSPs based on the Padova stellar models, calculated using the CMD 3.4 online tool ⁴ . We used PARSEC evolutionary tracks (Bressan et al. 2012; version 1.2S) for pre-main sequence to first thermal pulsation or carbon ignition, and COLIBRI models (Pastorelli et al. 2020) for thermal pulsation-asymptotic giant branch evolution. Circ*mstellar dust properties for M and C stars were from Groenewegen (2006), long-period variability along the red giant branch and asymptotic giant branch phases was taken from Trabucchi et al. (2019), and the initial mass function of Kroupa (2001) was used.

The goal is to model the broadband FUV–NUV-optical-IR SED of the M31 bulge. Because the broadband includes significant contributions from stars, gas and dust, a modeling program that includes these components is needed. Here we use the galaxy SED fitting program CIGALE (Boquien et al. 2019; Noll et al. 2009; Burgarella et al. 2005).

The fluxes for each annulus in the 17 filters listed in Table 1 were used as input into the CIGALE program for SED fitting. Various test cases were used to determine the best region shapes to use. For regions, we tested a series of boxes radially outwards from the nucleus, one large quarter circle region, multiple quarter circles, circular annuli, and elliptical annuli. The set of elliptical annuli was expected to best follow any spatial variation in the SED of the bulge of M31 because the annuli closely follow regions of constant surface brightness of M31's bulge (Leahy et al. 2021c). We confirmed that the CIGALE fits, as determined by χ², were as good or better for the elliptical annuli than other shapes of regions.

For SSPs, CIGALE uses the model of Bruzual & Charlot (2003) with Chabrier initial mass function. There are several options for star formation histories to use for combining SSPs of different ages. The nebular emission model has four parameters which we kept fixed at default values. These are ionization parameter, fraction of Lyman continuum photons escaping the galaxy, fraction of Lyman continuum photons absorbed by dust, and line width.

The dust attentuation model was a modified Calzetti et al. (2000) attenuaton law with eight parameters. These are: E(B − V)_lines, the color excess of the nebular lines; the reduction factor E(B − V)_factor to apply on E(B − V)_lines to compute E(B − V)_stars for the stellar continuum attenuation; λ_{uv bump}, the central wavelength of the UV bump; w_{uv bump}, the width of the UV bump; A_{uv bump}, the amplitude of the UV bump relative to that for the Milky Way; δ_powerlaw, the slope of the power law modifying the attenuation curve; ExtLaw_lines: the extinction law to use for attenuating the emission lines, either 1 for Milky Way, 2 for LMC or 3 for SMC; and R_V : the ratio of total to selective extinction.

CIGALE provides "Bayes" estimates and "Bayes" errors for each parameter requested by the user, summarized as follows (details given in Boquien et al. 2019). A large grid of user specified models is calculated (∼10⁵ to 10⁶) each with a different parameter set. For each model the χ² from the data versus model comparison and the likelihood () are calculated Then the Bayes estimate is the likelihood-weighted average and the Bayes error is the square root of the likelihood-weighted standard deviation for each parameter. For all of the results (parameter values) presented here from CIGALE, we use the Bayes estimate as the "best-fit" value and the Bayes error as the error in the "best-fit" parameter.

Initially we varied all parameters in a large number of CIGALE runs, each with ∼10⁵ parameter sets. For five of the parameters, the Bayes estimate was equal to the default initial value in CIGALE. These parameters and their Bayes values were: λ_{uv bump} (=217.5 nm); w_{uv bump} (=35.0 nm); A_{uv bump} (=1.0), ExtLaw_lines (=1); and R_V (=3.0). For subsequent CIGALE runs, these five parameters were fixed at these values.

Dust emission was modeled with the Dale et al. (2014) dust emission templates. This had two parameters frac_AGN the AGN fraction, and α_dust the slope of the dust emission. We tested varying frac_AGN but found that frac_AGN = 0 always worked best.

Tests were carried out to determine the best model star formation history of those available in CIGALE. We tested four different SFH models to find the best fit to the data, including rectangular periodic pulses (sfhperiodic), two decaying exponentials (sfh2exp), a delayed SFH with an exponential burst. The best model for the bulge of M31 was the delayed SFH with an exponential burst.

The delayed SFH had five parameters in CIGALE (Tau_main: the e-folding time of the main stellar population; age_main: the age of the main stellar population; Tau_burst: the e-folding time of the late starburst population; age_burst: the age of the late starburst population; and f_burst: the mass fraction of the late starburst population). The stellar emission model in CIGALE takes metallicity, Z, as an input parameter which is the same for all SSPs.

For the elliptical annuli, the fitting was done for each quarter of a given annulus and for the total flux from each annulus. The fit results for a given annulus between the different quadrants were consistent with eachother within errors in all case. Thus we present results only for fluxes from full annuli.

A set of parameters for each model component (star formation history, dust attenuation model, dust emission model) are input into CIGALE. Then the code evaluates a large number of models (∼10⁵ –10⁶) with one model for each combination of input parameters. The model with the smallest χ² is found and labeled as the best-fit model. CIGALE gives the SED for each region calculated with the parameters that gave the best fit for that region.

CIGALE also provides a Bayesian best estimate of parameters requested by the user, as noted above. In order for the Bayesian estimate to be reliable, the parameter spacing must be small enough that χ² does not change too much between adjacent parameter values, and it must be large enough that the largest χ² is larger than the smallest χ² by enough to give a large variation in . Thus, a number runs of CIGALE with different parameter spacings for each of the parameters were done. in order to obtain reliable Bayes estimates.

The CIGALE models produced good fits to the wide band FUV to IR SED of the 10 annuli and of the total bulge. An example of a best-fit model for the innermost annulus (Ann1) is shown in Figure 3. The data points for the five shortest wavelengths are from UVIT, the next five are from SDSS, the two points at 3.6μ and 4.5μ are from Spitzer/IRAC, and the five far-IR points are from Herschel/PACS and Herschel/SPIRE (Table 1).

The CIGALE models have a large number of input parameters. A few of the parameters were found to be consistent with constant values during test runs and were set to be fixed during the analysis runs of CIGALE, as noted in Section 3.2. The remaining parameters were free for our analysis. The analysis of the FUV to IR SED of the M31 bulge using the CIGALE models produced statistically acceptable models with χ² per degree of freedom of ≃1.

We use the Bayesian estimates from CIGALE as the best estimates of the parameters, noting that these are stable for a grid of input parameters with sufficiently many and closely-enough spaced parameters. Hereafter we label the Bayesian parameters as best-fit parameters. We divided the best-fit parameters into two categories: those without spatial dependence (independent of annulus number) and those that are spatially dependent. A parameter was considered not spatially dependent if the best-fit values for different annuli agreed with eachother to within the 2σ errors.

The parameters that do not have spatial dependence are listed in Table 2. The ratio of stellar extinction to line extinction and the dust emission power-law value are constant. The age and e-folding timescale for the main/old stellar population are constant, showing a rather old (12 Gyr) and rapid star formation episode (e-folding timescale of 450 Myr) for the main population of the bulge of M31. For the late star formation burst, the e-folding time is constant at ∼25 Myr. The metallicity, Z, is constant between different annuli and approximately solar metallicity. The errors on metallicity from CIGALE are not accurate because the allowed values of metallicity are chosen from a course (discrete) grid, making it not possible to compute fits with a closely-enough spaced grid of metallicity for the Bayesian estimate to be accurate. However, the course grid limits the metallicity to be between 0.008 and 0.05.

The best-fit spatially-variable parameters and errors are given in Table 3 for the 10 annuli and for the sum of the annuli (Total). R_eff is the value of R weighted over the area of the annulus.

Table 3.CIGALE Model Best-fit Parameters for the M31 Bulge FUV–NUV-optical-NIR-IR Data That Show a Radial Dependence

Region	Reff	Rout	χ2	E(B–V)lines	E(B–V)stars	δpowerlaw	error	error
				error	error	error	(1042 erg s−1)	(108 M⊙)
	('')	('')
Ann1	39.6	59.4	6.83	0.045	0.0188	−2.22	5.17	26.5
				0.012	0.0022	0.20	0.26	1.7
Ann2	72.3	83.9	10.74	0.065	0.0217	−2.16	2.92	14.7
				0.026	0.0014	0.20	0.15	1.1
Ann3	93.6	102.6	10.77	0.107	0.0293	−1.75	2.26	11.30
				0.034	0.0026	0.22	0.11	0.82
Ann4	110.9	118.8	7.25	0.114	0.0344	−1.50	1.938	9.63
				0.041	0.0030	0.21	0.097	0.71
Ann5	125.8	132.5	11.90	0.112	0.0334	−1.54	1.604	7.96
				0.038	0.0027	0.22	0.080	0.58
Ann6	139.1	145.4	10.99	0.114	0.0351	−1.44	1.491	7.35
				0.041	0.0036	0.22	0.075	0.54
Ann7	151.3	157.0	11.90	0.112	0.0409	−1.34	1.324	6.41
				0.035	0.0041	0.21	0.066	0.47
Ann8	162.4	167.8	11.77	0.119	0.0448	−1.21	1.202	5.88
				0.033	0.0049	0.19	0.060	0.43
Ann9	172.9	177.8	5.77	0.133	0.0515	−1.16	1.111	5.44
				0.028	0.0057	0.17	0.056	0.39
Ann10	182.7	187.6	8.19	0.135	0.0558	−1.12	1.052	5.14
				0.025	0.0058	0.15	0.053	0.37
Total		187.6	8.26	0.111	0.032	−1.59	20.2	101.3
				0.036	0.003	0.22	1.0	7.2

Note. a: The inner major axis for Total and Ann1 is 0'', for Ann2 to Ann10 the inner major axis is the outer major axis of the Annulus of one lower number.

Download table as: ASCIITypeset image

In summary (Table 2), the CIGALE analysis shows a dominant 12 Gyr old populations in the M31 bulge and a "600 Myr" young subdominant (f_burst ≃ 0.003) population. The masses of both populations and their extinctions vary with radius (Table 3).

As a first step, we verified the result from Leahy et al. (2021c) that the FUV–NUV SED of the bulge could not be fit with a single SSP. We found that the annuli could be fit satisfactorily with two SSPs. The best-fit two SSPs here have similar properties to those found by Leahy et al. (2021c): one SSP with large mass is an old population with low metallicity and the second SSP with small mass is young with higher metallicity.

The SSP models for the FUV–NUV 5 band data and the CIGALE models for the 17 band data yield very different model parameters. This is likely caused by the degeneracy inherent in fitting the five band FUV–NUV data with a model with more than five parameters. The degeneracy can be removed by including more data, thus we include here the optical and NIR data with the FUV–NUV data (total 12 bands) for the SSP modeling. We do not include the Herschel FIR data for which the dust and gas contributions, which are missing in the SSP models, are dominant (see Figure 3).

A single SSP model does not fit the FUV–NUV-optical-NIR SED of the bulge. The best-fit model reproduces the SDSS and UVIT F279N fluxes but is below the IRAC fluxes (about 30%) and well below the UVIT F148W, F169M, F172M and F219M fluxes.

The results of a two-SSP model applied to the FUV–NUV-optical-NIR SED of the bulge are given in Table 4. The 2-SSPs fit is considerably better than 1-SSP fit but is statistically rejected with high χ²(∼10⁴ to 10⁵). Figure 4 shows the best-fit 2-SSP model for Ann1 of the M31 bulge. A young SSP is required to fit the short wavelength UVIT data and an old SSP is required to fit the the SDSS data and the IRAC data.

Table 4.2-SSP Best-fit Parameters for the M31 Bulge FUV–NUV-optical-NIR Data

Region	Rout	χ2		age1		E(B–V)1		age2		E(B–V)2
	('')		(108 M⊙)	(Gyr)			(108 M⊙)	(Myr)
Ann1	59.4	1.4 × 105	22	4.0	0.30	0.17	9.1	560	−0.15	0.54
Ann2	83.9	1.2 × 105	20	4.0	0.30	0.23	1.6	380	0.30	0.33
Ann3	102.6	9.1 × 104	9.7	2.5	0.07	0.28	2.3	400	−0.12	0.59
Ann4	118.8	8.6 × 104	9.3	3.2	0.00	0.28	1.8	490	−0.15	0.51
Ann5	132.5	4.4 × 104	19	10	0.30	0.18	1.3	960	−0.42	0.25
Ann6	145.4	4.2 × 104	18	10	0.30	0.18	1.1	940	−0.34	0.21
Ann7	157.0	2.7 × 104	4.0	10	0.30	0.14	3.6	1100	−0.87	0.48
Ann8	167.6	3.6 × 104	14	10	0.30	0.18	0.91	970	−0.44	0.23
Ann9	177.8	3.0 × 104	9.3	10	0.30	0.17	1.8	1000	−0.73	0.42
Ann10	187.6	3.2 × 104	8.4	10	0.30	0.16	1.7	1100	−0.81	0.41
Total	187.6	3.7 × 105	290	10	0.30	0.20	8.6	490	0.30	0.18

Note. a: Subscripts 1 and 2 refer to the two SSPs: SSP₁ and SSP₂.

Download table as: ASCIITypeset image

We find that inclusion of optical and NIR data changes the best-fit models significantly. One difference from the fits to the UVIT only data is that the FUV–NUV-optical-NIR SED results in a significantly larger mass (factor ∼3) for the old population. Another significant difference is that the metallicity of the old SSP is found to be supersolar (+0.3), instead of the previous result of metal poor (−1.5 for the old SSP). The inclusion of optical and NIR data yields more reliable parameters than those from the NUV-FUV only fits.

We applied a three-SSP model to the FUV–NUV-optical-NIR SED of the bulge to see if this gave improved fits. The three-SSP model resulted in improved fits for Ann1, Ann2, Ann3, and Ann4 but not for the outer annuli, nor the total bulge. An example fit for Ann1 is given in Figure 5. For these inner annuli, the results are given in Table 5. The main differences are as follows. The extinction is lower for the three-SSP fits and closer to that from the CIGALE fits. For the three-SSP fits, the old population is now 10 Gyr age for the inner four annuli, the same as for the outer annuli from the two-SSP fits. The metallicity of the intermediate age population is now close to solar, and the age of the intermediate population is 600 Myr for all four inner annuli.

The CIGALE models include emission from stars, dust and gas, which make important contributions to the broadband SED. These three contributions can be seen for the best fit for Ann1 shown in Figure 3. For the fit shown in Figure 3, the contribution of gas is several percent for the FUV–NUV bands and ∼1% or less for the optical and NIR bands, and the contribution of dust is ∼1% for the NIR. Because the UVIT and SDSS photometry has small errors (≃0.05%–0.5%), these contributions are important to include in models.

Other advantages of the CIGALE models for fitting are that the dust attenuation law (see Section 3.2) has several adjustable parameters, and there are a number of SFH models to choose from (instead of the instantaneous bursts of the SSP models). Disadvantages of CIGALE include: stars, dust, and gas all have the same metallicity; all stars of different ages have a single metallicity values; and there is only a course grid of metallicity values to choose from.

Overall, the CIGALE models should be better than the SSP models, in particular for the SFH and dust and gas properties of the M31 bulge. Because the SFH is more realistic from CIGALE, the stellar masses from CIGALE should be a better measure of the true stellar masses. The disadvantages regarding metallicity mean that the metallicities from the SSP models are likely better.

For the two-SSP and three-SSP fits we found the three-SSP model gives better fits for Ann1 to Ann4, but not for Ann5 to Ann10, thus we use three-SSP results for Ann1 to Ann4 and two-SSP results for Ann5 to Ann10. A further justification of this is found in Dong et al. (2018), who find three stellar populations inside ∼100'' and 2 stellar populations outside, consistent with the SSP results here.

The extinctions for the SSP model fits are higher (∼0.15–0.45) compared to those for the CIGALE fits (∼0.11 for gas, ∼0.04 for stars). The SSP ages of the old populations (10 Gyr) are similar to the CIGALE ages (12 Gyr). The ages of the "600 Myr" population are very similar for SSP and CIGALE fits. The masses of the old population from the SSP fits are higher (by factor ∼2) than those from the CIGALE fits, likely because of the higher extinction of the SSP fits. The metallicities of old population and "600 Myr" populations for the SSP models are log(Z₁/Z_⊙) ≃ 0.3 and ∼0.0, compared to the single metallicity of log(Z₁/Z_⊙) ≃ 0 for the CIGALE fits.

In summary, the CIGALE and SSP fits agree on some results and disagree on others. The differences are mainly from two effects: the more flexible extinction law in CIGALE and the more flexible metallicities in the SSP models. From the CIGALE fits, the M31 bulge has a main stellar population which formed 12 Gyr ago, rather rapidly (with e-folding time of 450 Myr; Table 2), and a late burst of star formation with ages 500–700 Myr. From the SSP fits, the metallicity of the oldest population is significantly above solar and that of the "600 Myr" population is nearly solar. The SSP fits find that the inner bulge has an additional very young population with small mass and subsolar metallicity.

Many stellar population parameters from the CIGALE fits are found to be constant, independent of distance from M31's nucleus (center of the bulge). These are listed in Table 2. In particular, the age of the old population, the age of the "600 Myr" population (age_burst) and mass ratio of "600 Myr" to old population (f_burst) are independent of radius. Radially-variable properties are nebular (lines) extinction (Figure 6 left), stellar extinction (Figure 6 right), the dust attenuation law (Figure 7), and mass of the old population (Figure 8).

The radial dependence of the SSP parameters are compared to those from CIGALE. _Mold from the SSP fits (Figure 9 left panel ⁵ ) shows a similar decrease with radius as that from CIGALE but is higher by a factor of ∼2–3. The mass fraction of the "600 Myr" population from the SSP fits varies between ∼0.04 and ∼0.1 (Figure 9 right panel), higher than f_burst ≃ 0.003 from the CIGALE fits. The higher old population mass and higher "600 Myr" mass fraction are likely caused by the higher extinctions and higher old population metallicity for the SSP fits.

The age of the "600 Myr" population from the SSP fits (Figure 10 left panel) for the inner bulge (<120'') agrees with that from CIGALE, indicating that this parameter is measured reliably. For the outer bulge, the SSP fits yield a higher age than the CIGALE fits, which can likely be attributed to the higher extinction in the outer bulge for the SSP fits.

The metallicity from the CIGALE fits is solar and constant with radius. ⁶ The-SSP fits yield a nearly constant metallicity versus radius for the old population with value of [Fe/H] ≃ 0.3 (right panel of Figure 10). For the "600 Myr" population, [Fe/H]≃0 in the inner bulge then decreases to ≃−0.6 for the outer bulge. The youngest (25 Myr old) population exists only in the inner bulge and has subsolar metallicity.

Next we compare the results of this work with previous studies of the SFH and metallicity of the M31 bulge. Hammer et al. (2018) using hydrodynamical simulations in comparison with observations of M31, found that a 4:1 major merger 2–3 Gyr ago could explain many of the observed large-scale features of M31. The Pan-Andromeda Archaeological Survey of M31 was analyzed by McConnachie et al. (2018) for substructures in the halo. They showed that the most distinctive substructures were produced by a minimum of five different accretion events in the past 3–4 Gyr.

Saglia et al. (2018) derived SSP properties from Lick/IDS absorption line indices finding >80% of stellar ages were greater than 10 Gyr in the bulge and the stars were metal rich ([Z/H] ∼ 0.3 in the central 100''). Dong et al. (2018) use observations of the central 55 of M31 using HST ACS. They use CMD fitting of resolved stars to find that most stars (>70%) are older than 5 Gyr and have [Fe/H] ∼ 0.3. There is a second component of stars of age ∼1 Gyr over the whole bulge, and the central 130'' has a small fraction of stars that are <500 Myr old.

Previous works have found an old and metal-rich population in the bulge, although with less detail than Dong et al. (2018), including Dong et al. (2015), Saglia et al. (2010), and Olsen et al. (2006). Evidence of younger stars in the bulge has been given previously. E.g. Dong et al. (2015) find a 300 Myr to 1 Gyr population with 1% of the mass. Dong et al. (2018) finds CMD fitting has better age resolution than other methods, which allows identification of separate 1 Gyr and <500 Myr populations.

The SSP fits presented here have an old (∼10 Gyr) metal-rich SSP consistent the with old metal-rich population previously identified. The CIGALE models find the age to be 12 Gyr with solar metallicity, where the solar metallicity is probably an artifact of having a single metallicity for the populations in CIGALE. Figure 21 and Tables 1 and 2 of Saglia et al. (2018) show an age of the old population of ≃11 to 13 Gyr nearly independent of radius (along the bar minor axis), in agreement with the CIGALE results. Their composite stellar models (their Table 4) interpret the age and metallicity as mixtures of four populations, a classical bulge of age 13 Gyr with [Z/H] = 0.35, a bar with age 12 Gyr and [Z/H] = 0, a boxy-peanut with age 13 Gyr and [Z/H] ≃ −0.1, and a disk with age 7 Gyr and [Z/H] = 0.01. Dong et al. (2018) find age of the old population to be >5 Gyr and [Fe/H] ∼ 0.3 at all radii. The 12 Gyr age of the old population from the CIGALE fits and the metallicity of the old population, [Fe/H] = 0.3 from the SSP fits are consistent with both studies.

For the "600 Myr" population, Figure 21 of Dong et al. (2018) shows an age of 0.9–1.5 Gyr, increasing with radius, and metallicity of [Fe/H] of 0.4–0.1, decreasing with radius. ⁷ A "600 Myr" age constant with radius was found by Dong et al. (2015) by modeling photometry of the bulge, in agreement with the current CIGALE and SSP results. Similarly, that work found a metallicity for the "600 Myr" population ([Fe/H] ≃ 0) in agreement with our work. The CMD analysis (Dong et al. 2018) yields higher metallicity and age than the photometry analysis (this work and Dong et al. 2015). Overall, we conclude the intermediate age population has age "600 Myr" to 1 Gyr and metallicity [Fe/H] ∼ 0.0–0.3, with refinement depending on development of better models.

The mass fraction of the "600 Myr" population was found to be constant in the CIGALE model at ≃0.003 and in the SSP model at ∼0.05, whereas Dong et al. (2018; their Figure 27) find values of 0.08–0.03. The SSP result agrees with that from Dong et al. (2018), indicating that the restrictions of CIGALE (same extinction and metallicity of both populations, and no third population) are skewing the CIGALE values.

The mass fraction and the age of the youngest population from Dong et al. (2018) are larger than that deduced from the SSP model, however their data did not include any FUV photometry so is less sensitive to the youngest populations. Because a lower mass is required for a young (brighter) population, their results are approximately consistent with the SSP results. The inclusion of FUV data here likely means the SSP result is a better measure of the age and mass of the youngest population.

The "600 Myr" and young populations could be the result of the multiple merger events over the past 3–4 Gyr, previously identified. Such mergers could yield the gas that reached the bulge and formed stars there ∼600 Myr and ∼25 Myr ago. The "600 Myr" and young stars would be the same metallicities, respectively, as the infalling gas, so that the gas was solar abundance (600 Myr) or subsolar (25 Myr). Alternately, the gas was not directly from the merging mass but from within M31 where its infall was triggered by the merger. For that case too, the infalling gas would have to have solar abundance (600 Myr) or subsolar (25 Myr) abundance prior to its infall and star formation in the bulge.