The Star Formation History in the M31 Bulge

The Star Formation History in the M31 Bulge

Hui Dong, Knut Olsen, Tod Lauer, Abhijit Saha, Zhiyuan Li, Ruben García-Benito, Rainer Schödel Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomá S/N, E-18008 Granada, Spain National Optical Astronomy Observatory, Tucson, AZ, 85719, USA School of Astronomy and Space Science, Nanjing University, Nanjing, 210093, China Key Laboratory of Modern Astronomy and Astrophysics at Nanjing University, Ministry of Education, Nanjing 210093, China E-mail:

We present the study of stellar populations in the central 5.5′ (1.2 kpc) of the M31 bulge by using the optical color magnitude diagram derived from HST ACS WFC/HRC observations. In order to enhance image quality and then obtain deeper photometry, we construct Nyquist-sampled images and use a deconvolution method to detect sources and measure their photometry. We demonstrate that our method performs better than DOLPHOT in the extremely crowded region. The resolved stars in the M31 bulge have been divided into nine annuli and the color magnitude diagram fitting is performed for each of them. We confirm that the majority of stars (70%) in the M31 bulge are indeed very old ( 5 Gyr) and metal-rich ([Fe/H]0.3). At later times, the star formation rate decreased and then experienced a significant rise around 1 Gyr ago, which pervaded the entire M31 bulge. After that, stars formed at less than 500 Myr ago in the central 130″. Through simulation, we find that these intermediate-age stars cannot be the artifacts introduced by the blending effect. Our results suggest that although the majority of the M31 bulge are very old, the secular evolutionary process still continuously builds up the M31 bulge slowly. We compare our star formation history with an older analysis derived from the spectral energy distribution fitting, which suggests that the latter one is still a reasonable tool for the study of stellar populations in remote galaxies.

Keywords: galaxies: bulges Galaxies, galaxies: evolution Galaxies, galaxies: individual (M31) Galaxies

1 Introduction

The surface brightness of bulges in spiral galaxies is smooth and their morphologies are similar to elliptical galaxies. Therefore, at first, bulges were thought to form through the same mechanism as elliptical galaxies, such as rapid collapse in the early Universe (Sandage, 1990) or major mergers (Toomre, 1977), which remove molecular clouds, stop star formation and increase velocity dispersion. In both cases, old stars (at least several Gyr old) are expected. In the former case, the stellar metallicity should be very low, while the metallicity in the latter cases depends on when the major mergers happened. On the other hand, young stars have also been found in bulges and are explained by another mechanism, secular evolutionary process, which could also play a significant role in shaping bulges; Materials in the disks spiral into inner regions of galaxies, assisted by the torque of galactic bars, trigger star formation and produce young stars (see the excellent review given by Kormendy & Kennicutt, 2004). Key to disentangle these mechanisms is to study the spatial distribution of the age and metallicity of stars in spiral bulges.

The Andromeda galaxy (M31) is the closest external grand design spiral galaxy (780 kpc, McConnachie et al. 2005; de Grijs & Bono 2014) and has a distinct bulge (Beaton et al., 2007), as well as a super-massive black hole (M31*, Dressler & Richstone, 1988; Kormendy, 1988; Crane et al., 1992; Garcia et al., 2010; Li et al., 2011). The Sérsic index of the M31 bulge at the band is 2.2 (Courteau et al., 2011). Therefore, the M31 bulge is classified as a classical bulge according to the criterium given in Fisher & Drory (2010). On the other hand, Beaton et al. (2007) identified a thin bar in the M31 bulge (see also Blaña Díaz et al. 2017; Opitsch et al. 2018), which is the major characterization of pseudobulges (Kormendy & Kennicutt, 2004). Compared to our own Milky Way, M31 is a better lab to study the formation and evolution history of bulges for two reasons: First, the M31 bulge suffers from small foreground extinction from the Milky Way (0.166 mag, Schlafly & Finkbeiner 2011), because of its Galactic latitude (-21.57 degree), and local extinction from a small amount of molecular clouds in the M31 bulge (Melchior & Combes, 2013; Dong et al., 2016). Therefore, unlike the Galactic bulge, the systematic uncertainty introduced by extinction can be ignored. Second, the M31 bulge is distant enough that we do not need to consider the differential distances of individual stars, while the M31 bulge is still close enough that we can resolve individual bright stars.

The biggest challenge to characterize the stellar populations in the M31 bulge is the crowding due to high surface brightness (Rich et al., 1993; Depoy et al., 1993; Rich & Mighell, 1995; Davidge, 2001; Stephens et al., 2003). Stephens et al. (2003) analyzed the observation taken by the Hubble Space Telescope (HST) of the M31 bulge with an angular resolution of 0.2″(0.8 pc). They found that their stars were significantly dimmer ( 1 mag) than the ones reported in Rich & Mighell (1995), which used the HST dataset with the original aberrated optics, and in Davidge (2001), which used the 3.6 m Canada-France-Hawaii Telescope with a seeing of 0.35″ (1.3 pc). This magnitude difference was apparently caused by the blending of multiple sources in the observation with poor angular resolutions.

In order to overcome the crowding problem, Olsen et al. (2006) (see also Davidge et al. 2005) analyzed adaptive optic (AO) images of the M31 bulge taken by the Gemini/NIRI camera in the near-IR (NIR) band, with an angular resolution of 0.1″ (0.4 pc). Although their observation was still not deep enough to reach red clump (RC) stars, they detected a significant amount of red giant branch (RGB) stars and asymptotic giant branch (AGB) stars. Through fitting of the NIR color magnitude diagram (CMD), they found that stars with 6 Gyr and solar metallicity contribute significantly to the M31 bulge (60%). Saglia et al. (2010) analyzed optical low resolution spectra taken by Hobby-Eberly Telescope and derived similar results, while they suggested an even older age (12 Gyr). They also found that the stellar metallicity decreased from 3 in the nucleus to 1 at 2′ (454 pc) away from the centre. Dong et al. (2015) performed the spectral energy distribution (SED) fitting of the photometric data derived from HST images at ten bands from ultraviolet (UV) to NIR and obtained similar results. Meanwhile, no stellar population younger than 10 Myr seems to exist in the M31 bulge (Rosenfield et al., 2012).

On the other hand,  Bender et al. (2005) reported that the HST/UVIS spectra of the UV bright spot identified by King et al. (1995) in the M31 nuclear region is similar to the one of A0-type stars. They agreed with the scenario proposed by Lauer et al. (1998) and Brown et al. (1998) that the UV bright spot could represent a 200 Myr old stellar cluster. Lauer et al. (2012) further resolved a handful of UV-bright stars within the central 0.2″ (0.8 pc), which were suggested to be young early-type main-sequence (MS) stars instead of old post early AGB (PEAGB) stars. Saglia et al. (2010) also found that the age of a single stellar population fitting the observed Lick indices decreased to 8 Gyr, inside the central 2″ (8 pc), which was much younger than that beyond the nuclear region (12 Gyr); An extra young stellar population with age 600 Myr and mass fraction lower than 10% could remove this age difference of old stellar population inside and beyond 2″. They suggested that a burst of star formation could have occurred within the nucleus 100 My ago.

That raised the question whether this star formation activity and the corresponding young MS stars are constrained to only the central few arc seconds or not. Dong et al. (2015) claimed that intermediate-age stellar population (300 Myr to 1 Gyr) could pervade the M31 bulge, while contributing only 1% of the total stellar mass, which is also supported by Conroy & van Dokkum (2016). However, this question has not been completely settled, because previous studies with integrated spectra or photometry could be dominated by a few bright stars, such as AGB stars in the NIR band or numerous old evolved low mass stars, such as RGB stars in the optical band. Even in the UV, the integrated intensity could still suffer from contamination of the numerous extreme horizontal branch stars (Dong et al., 2015, and references therein). Meanwhile, limited by the total number of photometric data points, only one or two starburst activities could be used to fit the whole star formation history (SFH) in the M31 bulge (Saglia et al., 2010; Dong et al., 2015). Therefore, we can not have a clean scrutiny of the SFH in the M31 bulge. On the other hand, CMD fitting technique uses all individual resolved stars and the result is not dominated by a few bright stars. Although the AO NIR images taken by ground-based 8m telescopes could indeed achieve diffraction limit, such as Olsen et al. (2006), the limited NIR color range (e.g., -) hampers our attempt to disentangle stellar populations with different ages and metallicities. Instead, because MS stars are brighter in the UV and optical bands, the corresponding CMD could be better to identify the intermediate-age stellar populations.

In this paper, we will utilize new deep UV and optical images taken by HST to perform CMD fitting and further study the SFH of the M31 bulge. The angular resolution of HST in the UV and optical bands is 0.05″, which is better than previous AO observation in the NIR band. On the other hand, HST has a more stable PSF, a lower background and a larger field-of-view (FoV) than the ground-based AO observation. All of these properties above are important for us to 1) detect dim stars, 2) reduce the photometric uncertainties in the crowded environments, which are critical for the CMD fitting, and 3) study spatial variations of different stellar populations.

The structure of this paper is as follows. In § 2, we present our observation. Then, we describe our methods to produce the Nyquist-sampled images, construct catalogues from detected sources in deconvolved images, perform artificial star tests, compare with the source catalogs from the pipeline of the Panchromatic Hubble Andromeda Treasury (PHAT) survey and discuss the foreground and background contamination in our catalog in § 3. In § 4, we perform the fitting of the CMD to disentangle stellar populations with different ages and metallicities in the M31 bulge and examine the blending effect on our CMD fitting. In § 5, we discuss the impact of the results on our understanding the assembly history of the M31 bulge. We will summarize our results in § 6.

2 Observations

We utilize the dataset taken by the HST/ACS WFC and HRC cameras. The FoV of these two cameras are 202″202″ and 29″26″, respectively. The former observation is from the PHAT survey (Program GO-12055, Dalcanton et al., 2012), while the latter one is from Program GO-10571 (Lauer et al., 2012). The detailed description of these two dataset are given in Dalcanton et al. (2012)Williams et al. (2014) and Lauer et al. (2012). Here, we only briefly summarize the key points, related to this work.

2.1 PHAT survey

The PHAT survey used the F475W and F814W filters in the optical band to satisfy the high throughput and large wavelength leverage simultaneously. At these two bands, each pointing has five and four dithered exposures, respectively. The basic information of the filters and exposure times of dithered images is listed in Table 1. The short ‘guard’ exposures (10 second in F475W and 15 second in F814W) were used to recover the photometry of bright stars, which saturated in the long exposures. Because the plate scale of the WFC camera (0.05″ pixel) is similar to the full width half maxim (FWHM) of the point spread functions (PSFs) at these two bands, the PSF is undersampled. In order to overcome this problem, the survey used half pixel dithering in four and two dithered images with the longest exposure times at these two bands, respectively. In §3.1, only these exposures are used to produce the Nyquist-sampled summed image. As a result, without the ‘guard’ exposures, we cannot recover the photometry of stars with F475W 18.7 mag and F814W 18.6 mag.

We used the data of Brick 1, one of the 23 Bricks in the PHAT survey. Each Brick consisted of 36 gird of pointings; each pointing was dubbed ‘Field’ (Dalcanton et al., 2012). We only concentrated on the half of Brick 1 (33 grid of pointings), which includes M31* and the M31 nuclear bulge; unfortunately, the other half, did not operate with the half-pixel dithering strategy, which was required to produce the Nyquist-sampled images below. The region used in this study is illustrated in Fig. 1a and the total sky area is 61 arcmin (870 pc). Besides the M31 bulge, the celestial northwest part of our FoV (see Fig. 1a) covers a small portion of the 5 kpc star formation ring (Lewis et al., 2015).

2.2 Hst ACS/HRC observation

Besides the PHAT dataset, the M31 nucleus was observed by the HRC camera at the F330W and F435W bands in the shorter wavelengths, with a sightly higher angular resolution (PSF FWHM 0.04″). Half-pixel dithering pattern was used in the observation. The pixel scale of the HRC camera (0.028″0.025″) is finer than that of the WFC camera by a factor of two. As a result, compared to the WFC observation, the PSF in the HRC images is better-sampled. Meanwhile, the two bands of the HRC observation in the near-UV provide us with a better chance to discover the intermediate-age stellar population, which is supposed to be brighter than the dominant old metal-rich population at this wavelength range. Therefore, the HRC observation is a good complement of the PHAT survey. Table 1 gives the central wavelengths and widths of these two filters, as well as the exposure times. The saturation magnitudes of these two bands are 13.6 and 14.8 mag. The final HRC F330W Nyquist-sampled image derived in §3.1 is shown in Fig. 2.

2.3 Spitzer/Irac 3.6 m observation

Besides the HST/ACS observations, we used the Spitzer/IRAC 3.6 m observation of the M31 bulge (Barmby et al., 2006) to constrain the spatial distribution of the stellar mass in the M31 bulge. As mentioned in §4, in the inner region, limited by our detection limit, the CMD cannot reach RGB stars older than 2 Gyr; As a result, we cannot constrain the stellar mass of the old stellar population, which is the major contribution of the total stellar mass. On the other hand, at 3.6 m, the stellar light is less affected by the foreground extinction and the corresponding mass-to-light ratio is insensitive to the age and metallicity of underlying stellar population (McGaugh & Schombert, 2014; Meidt et al., 2014)McGaugh & Schombert (2014) and Meidt et al. (2014) give 0.47 and 0.6 M/L, with Kroupa and Chabrier initial mass function (IMF) (Kroupa, 2001; Chabrier, 2003), respectively. In this paper, we choose the value given in McGaugh & Schombert (2014) to convert the 3.6 m intensity map into the stellar mass map, because we use Kroupa IMF in the CMD fitting below.

3 Data Reduction

3.1 Nyquist-sampled images

Lauer et al. (2012) have analyzed the same HRC dataset with a Fourier method, which was originally presented in Lauer (1999). We employed the same technique to reduce the WFC observation and produce Nyquist-sampled images for individual pointings. Unlike the widely used ‘AstroDrizzle’, the technique of Lauer (1999) does not cause any smoothing and degradation of the PSF. The detailed description of the procedure is given in Lauer (1999) and Lauer et al. (2012). We simply describe some key points below.

We constructed Nyquist-sampled images from calibrated dithered exposures kindly provided by the PHAT pipeline team. These exposures were bad pixels removed, bias and dark subtracted, as well as flat-fielding corrected. Charge transfer efficiency which the WFC camera suffers was amended. All the dithered exposures of the same filter and pointing were aligned. We first rescaled the calibrated images at the F475W and F814W bands in units of counts to have 360 and 800 seconds respectively. They were then interlaced into a roughly Nyquist-sampled image with a square-spiral pattern. The pixels with abnormal count rates were selected out and replaced by the interpolation from nearby pixels. Second, we used the Fourier algorithm of Lauer (1999) to combine the repaired dithered exposures into a well-sampled Nyquist images. This algorithm sums the Fourier transforms of the dithered exposures, suppresses the high-frequency power and transforms the result back into imaging space, but with half of the original WFC pixel size, i.e. 0.025″/pixel. Third, using the calibration files from the STScI, we corrected for the distortion and multiplied the pixel area map to get the final Nyquist-sampled images at different pointings and filters.

For each Nyquist-sampled image, we also produced a sky background and derived a noise value for source detection below. Individual Nyquist-sampled images were smoothed with a median box, the size of which is 250250 pixel (6.25″6.25″) and 100100 pixel (1.14″1.14″) for the WFC and HRC observations, respectively, to be their sky background images. The histogram of the pixel values of individual Nyquist-sampled images includes a Gaussian function plus a tail toward the bright end; The former one represents the fluctuation introduced by the bias and sky background and the latter one is contributed by bright stars. We derived the standard deviation of the Gaussian function and used it to be the noise of each Nyquist-sampled image.

Then, at each band, we merged these Nyquist-sampled images together to the final mosaic. We employed the method detailed in Dong et al. (2011) to derive the relative astrometry and bias offset among different pointings, by using their overlap regions. Then, we used bright isolated stars to align our images to the absolute astrometry system defined by the 2MASS catalog (Skrutskie et al., 2006). The F814W mosaic of our data is shown in Fig. 1b with the boundaries and labels of the ‘Fields’.

3.2 Deconvolution photometry

We first empirically constructed the PSFs from bright isolated sources. Considering the spatial variation of PSFs throughout the large FoV of the WFC camera, we derived the PSFs for the four quadrants of the two chips, respectively. In the M31 bulge, limited by the crowding, very few isolated stars in individual quadrants of one pointing are bright enough to produce the PSFs. Therefore, we aggregated isolated bright stars (20, 19.5) in the same quadrant of the nine pointings, considering that these images were taken during the same period (around half month) and the HST PSF was very stable. We removed faint point sources near these bright stars, normalized and median-averaged their images to obtain the final PSFs at different quadrants in the two chips, respectively (i.e. 42=8 PSFs for each band). For the HRC observation, considering its small FoV, only one PSF was produced for the F330W and F435W observations, respectively. Since numerous RGB stars are dimmer in these short wavelengths, the surface number density is low and there are still many isolated bright sources which could be used to produce the PSFs.

Then we used the empirical PSF derived above to deconvolve the Nyquist-sampled images using the Lucy-Richardson algorithm. This algorithm worked iteratively to reduce the wing of the PSF and enhance its core, simultaneously. Therefore, in the deconvolved images, individual stars became delta functions. We found that 320, 160 and 640 iterations were enough for the WFC, HRC F330W and F435W Nyquist-sampled images, respectively.

We detected sources and extracted their photometry from the deconvolved images. First, we selected the pixels in the deconvolved image with values larger than the background image by three times of the noise (the background images and the noises were produced in §3.1). We refer to these pixels as ‘source’ pixels. We sorted them according to their intensities. Second, starting with the brightest one, we considered that the adjacent 33 pixels belong to the same source and then produced a source list. Third, we merged the source catalog of the two bands (hereafter, ‘two bands’ refer to F475W and F814W for the WFC camera or F330W and F435W for the HRC camera) into a master catalog. Two sources at these two bands were considered as the same one, if their distance was less than 2 pixels, i.e. 0.05″(0.025″), in the WFC (HRC) observation. In this master catalog, we also kept the sources detected in only one band. Fourth, the intensity of each source in the master catalog at each band was the difference of the sum of values of its ‘source’ pixels in the Nyquist-sampled image and the background image. Fifth, we calculated the signal-to-noise ratio of each source as the ratio of its intensity and the square root of the sum of values of its ‘source’ pixels in the Nyquist-sampled image (which is equal to intensity plus background). Sixth, we translated the intensity into magnitude with the zero magnitudes: 22.909, 25.188, 26.157 and 25.523 mags for the F330W, F435W, F475W and F814W bands, respectively111 Seventh, for the WFC observation, we compared our source list of each position with that from the PHAT pipeline and derived the mean difference is only 0.05 mag. We subtracted this value from the source catalog of each pointing. In our final catalog, we only included the sources with signal-to-noise (S/N) ratio 3 in both two filters of the WFC and HRC observations. Meanwhile, because we did not use the guard exposures in the WFC observations when producing the Nyquist-sampled images, we could not recover the photometry for the bright sources. We identified a saturated source, if any of its adjacent 1111 pixels had count rate 109 (95) electrons/s at F475W(F814W) bands, according to the ACS Instrument handbook. These stars were removed away from our discussion below. Eighth, for the WFC observations we merged the source catalogs of the nine pointings into the final source catalog for the CMD fitting in §4.

3.3 Artificial star test

We generated artificial stars with -5F475W-F814W10 and 13F814W29 (-4F330W-F435W4 and 18F435W28), added them back into the original WFC (HRC) images and re-performed the source detection. Through comparing the input and output magnitudes, we determined the recovered fraction (i.e., the fraction of input artificial stars, which were detected) and uncertainties as a function of magnitude, color and local surface brightness. The dimmest magnitudes of the artificial stars were chosen to be one magnitude lower than the cut-off magnitudes of our source catalogs to examine the potential Eddington-bias.

We used the Padova stellar evolutionary tracks (version: PARSEC v1.2S +COLIBRI PR16, Marigo et al., 2013, 2017; Rosenfield et al., 2016) to generate the artificial catalog, which consists of two groups, both of which have similar numbers of stars: 1) a stellar population with continuous star formation, as well as various metallicities (from 1./200 to 2.5 ) and ages (from 4 Myr to 13 Gyr). This component reflects the color and magnitude distributions of all potential stars. 2) a stellar population with old age ( 5 Gyr) and high metallicity (). These red and dim stars dominate the M31 bulge and are near or below the detection limit of our observations (see below). As a result, only very few of these input stars in the previous group have been detected and the derived recovered fraction at this color and magnitude ranges in the CMD suffers from large statistic uncertainty. Therefore, we added a new old and metal-rich component in the artificial catalog to enhance the number of recovered stars near the detection limit.

For the WFC dataset, we manufactured ten observations at each pointing and band for the two groups above. In each manufactured observation, we injected 1 artificial sources with the PSF constructed in §3.2 into the original Nyquist-sampled image of the WFC observation, which was 5% of the stars found in each pointing. For the HRC observation, we produced 10 artificial observations, each including 350 artificial sources. We ran the deconvolution routine and performed the source detection. We recorded the positions, magnitudes and the galactocentric radii () of individual output artificial stars.

Besides the input magnitude, the recovered fraction and photometric uncertainties are determined by local surface brightness too. In the M31 bulge, the surface brightness is axisymmetric with respect to the M31 nucleus (Beaton et al., 2007; Dong et al., 2014, 2015) and changes quickly along . As a result, due to the large FoV, even in the same WFC pointing, the recovered fraction and photometric uncertainties could change at different locations of the camera. Thus, instead, we equally divided the artificial stars in all the pointings, according to their , into nine annuli. The inner and outer of these nine annuli are given in Table 2. When increases, the stellar number surface density decreases and our source catalog is less affected by the blending, especially for faint sources. Thus, the recovered fraction increases. Fig. 3 shows the recovered fraction as a function of the input magnitude at the nine different radii. The 50% detection limit of the F475W and F814W bands decreases from 26.3 and 24.4 mag for the outermost ring by 2 mag to 24.6 and 22.3 mag in the innermost ring. For the HRC observation, due to the small FoV, we did not divide it into finer parts. The 50% detection limit of the F330W and F435W bands are 24.6 and 23.9 mag, respectively. The recovered fractions of these two bands are also given in Fig. 3.

Fig. 4 shows the difference between the output and input magnitudes of the artificial star tests for the innermost and outermost radii of the WFC observation, as well as the HRC observation, as a function of the input magnitudes. The deviation for the bright stars is small. Instead, it increases towards faint magnitudes. Due to blending, the output magnitude is averagely smaller than the input one. For bright stars, the deviation is larger at the F814W band than the F475W, F435W and F330W bands, indicating the crowding issue is more significant at the long wavelength.

3.4 Comparison with the DOLPHOT photometry

We compare the performance of our deconvolution method and the DOLPHOT algorithm used by the PHAT pipeline (Dalcanton et al., 2012) in the extremely crowded environment, like the M31 bulge, through the real observation and artificial star clusters.

3.4.1 Real observation

Figs. 5 and 6 are the comparisons of the CMD and luminosity function (LF) between the DOLPHOT ‘st’ catalogs222The PHAT pipeline provides two catalogs for each ‘Field’: ‘gst’ and ‘st’. The former uses stricter quality control. As a result, the incompleteness is high in the crowded regions. Therefore, we used the ‘st’ catalog for comparison. and our ones at two ‘Field’s: Field 11 and Field 18. These two pointings have the highest and lowest number surface densities in our study. As a result, the detection limit at the F814W band in Field 18 is 1 mag below that in Field 11. The CMD and LF derived from the DOLPHOT and deconvolution method are overall very similar. However, compared to the DOLPHOT catalogs, our deconvolution method can detect more dim and red stars in both fields; for example, our deconvolution method detects roughly four times more stars below the detection limit of the F814W band than DOLPHOT. Because we required the sources to have S/N3 at both F475W and F814W bands, our LF has a clear cut at the dim sides of these two bands.

3.4.2 Artificial star clusters

We inserted 28 artificial stellar clusters into original observation of Brick 17 and Field 06. The salient parameters of these star clusters are listed in Table 3. Compared to the real observation, the advantage of this strategy is that we know the exact input magnitudes and positions of artificial stars. These clusters span a large range in the total mass and compactness, from the very massive and compact ones, such as, ID 15 and ID 19, to small and loose ones, for example, ID 8. Thus, we could compare the performance of DOLPHOT and deconvolution photometry at different environments.

We used both our deconvolution method and DOLPHOT algorithm with the parameters given in Dalcanton et al. (2012) to extract the photometry. Fig. 7 shows the WFC F814W images, the DOLPHOT and the deconvolved CMDs in the central 2″, with the field sources removed, of four representative star clusters. These clusters were sorted according to their mass densities. ID 19 is the most compact artificial star cluster. However, the large foreground extinction attenuates the brightness of individual stars. Therefore, most of the sources detected by the DOLPHOT and deconvolved method should be the blending of several faint sources and become brighter than the input magnitudes. On the other hand, ID 1, has much less attenuation. Thus, we can see an apparent compact star cluster in the F814W image. Interestingly, the DOLPHOT CMD shifts to the brighter side of the input CMD, while the deconvolved and input CMDs match better. The similar phenomena is found in the ID 16, which is less compact than ID 1, by a factor of 3. In conflict, for ID 6, the mass density of which is smaller than that of ID 1 by an order, both the CMDs derived by the DOLPHOT and deconvolution method are similar to the input one.

We further compared the input and output magnitudes and colors of individual stars to understand the difference between the catalogs derived by the DOLPHOT and deconvolution methods in Fig. 8. We chose two artificial star clusters, ID 1 and ID 6. The output magnitudes of dim stars are smaller than the input ones, due to the blending. Unlike the magnitude, the input and output colors are very similar (0.3 mag). That’s because the blending effect at both the F475W and F814W bands cancel each other. In ID 1, the deconvolved method performed better in the photometry than the DOLPHOT method, in terms of the smaller difference between input and output magnitudes. For example, the mean and standard deviation of the F814W magnitude are 0.90.93 mag for DOLPHOT and 0.570.83 mag for deconvolution. This explains why the output CMD from the DOLPHOT algorithm is brighter than that from the deconvolution method.

Fig. 9 shows the surface brightness (mag arcsec) at the F814W band as a decreasing function of (in units of arc second) in our FoV. The region inside the central 230″ is brighter than ID 1 (18 mag arcsec). Therefore, the deconvolution method could provide us with a better catalog to study of stellar population in the M31 bulge.

3.5 Foreground and background contamination

Because the galactic latitude of M31 is low, our detected sources should include foreground Galactic dwarf stars, which have 1F475W-F814W4. However, they account for a very tiny fraction of our observed stars. For example, we used the ‘Trilegal’ code (Girardi et al., 2012) to simulate the stellar population of the Milky Way in the FoV of our WFC observation. The code outputs 117 stars and only 113 of them are brighter than the 50% completeness limit at both the F475W and F814W bands of the outermost #9 annulus. The HRC observation should includes even much less contribution from the foreground Galactic dwarf stars, because of the smaller FoV and the faintness of dwarf stars at the UV bands.

Dalcanton et al. (2009) suggested the number density of background galaxies in the optical bands to be 60  (see also Williams et al., 2014). Therefore, we only expect less than one background galaxy in our FoV.

4 Analysis and Results

In §4.1, we first qualitatively analyzed different components in the CMD derived from our deconvolution catalog. Then we quantitatively fit the CMD to get the SFH in the M31 bulge in §4.2. After that, in §4.3, we examine the blending effect on our result based on the CMD fitting.

4.1 Cmd

The Hess representations of the CMD of our ACS observations are shown in Fig. 10 for both the WFC (top panels) and HRC (low panels) observations. We used 0.02 mag and 0.05 mag as the bin sizes of the color and magnitude terms in the CMD of the WFC observation, due to its large FoV, while 0.2 mag and 0.5 mag bin sizes are selected for the HRC observation. In the WFC observation, most stars concentrate in the regions enclosed in the red parallelogram, which are redder than most of the metal-poor isochrones and instead should be RGB stars of the metal-rich populations. On the other hand, the blue stars embedded by the blue parallelogram could be MS stars, and/or post-AGB stars (Rosenfield et al., 2012). Like Monachesi et al. (2011), we called this region in the CMD ‘blue plume’ (BP). In the HRC observation, most of the stars fall into the regions of the CMD close to or below the 50% detection limit. On the other hand, the stars brighter than F435W=24 mag are divided in two branches: one with F330W-F435W-1 and the other one with F330W-F435W1.

We further present the evolution of the CMD along for the WFC observation. Fig. 11 shows the Hess representations of the CMD at the nine annuli defined in Table 2. In each annulus, only stars brighter than 50% completeness limit are kept. Due to the crowding, the recovered fraction in the innermost annulus is lowest. As a result, from the CMD, it is impossible to constrain the star formation rate (SFR) for the stellar population with super-solar metallicity with age older than 2 Gyr. On the other hand, the RGB stars with super-solar metallicity and 10 Gyr old are still brighter than 50% completeness limit at both bands in annulus #9. The BP stars with F475W-F814W can be found in all the nine annuli, but mainly concentrate in annuli, #8 and #9 at F814W23-24 mag.

In order to further study the nature of stars embedded within the blue parallelogram in the top left panel of Fig. 10, Fig. 12 shows their spatial distribution. The majority of these stars concentrate in the M31 nucleus. Then, the surface number density decreases till 230″ away from M31*, then increases again, caused by the enhancement of stars in the northwest of our FoV. Fig. 13 shows the CMD at three different regions, divided by the two cyan circles with of 80″ and 230″. In the lower left panel of Fig. 13, the stars beyond 230″ concentrate near - and 22.5. These stars should be MS stars in the 5 kpc star formation ring. Therefore, #8 and #9 annuli given in Table 2 include not only the bulge stars, but significant disk contribution. On the other hand, the stars within the BP in the central 230″ distribute more uniformly in the CMD and do not seem to come from the young stellar population.

Fig. 14 shows the differences of the CMDs between the six inner annuli and #7 annulus to further demonstrate the evolution of the CMD in the M31 bulge. Unlike #8 and #9 annuli, these inner annuli should still be dominated by stars in the bulge. Before subtracting the CMD of #7 annulus, we first corrected the difference of the recovered fraction among annuli and then scaled the CMD of #7 annulus to the same total number of stars within the 50% completeness limit of the individual inner annuli. Compared to the #7 annulus, the innermost annulus has extra stars at F814W22 mag and F475W-F814W1.5, as marked by the cyan box. This region is redder than the stars in the BP, but can still be explained by the stellar population with solar metallicity and 200 Myr old. These stars become dimmer with the increase of from #1 to #4 annuli.

4.2 Fitting the CMD

We used the source catalogs given in §3.2 and the artificial star tests performed in §3.3 to constrain the stellar populations in the M31 bulge. Specifically, we utilized the CMD fitting method developed by Olsen et al. (2006) to compare the observed data with those predicted by the stellar evolutionary isochrones with different ages and metallicities. This method first uses the stellar evolutionary tracks to generate CMDs, giving a set of parameters (e.g. age, metallicity, distance, reddening, and initial mass function). The synthetic CMD is convolved with the observed photometric uncertainty to simulate the real observation. The likelihoods of the observed CMD as a representation of a set of starburst stellar populations are calculated and are used to determine the best solution.

The Padova stellar isochrones were utilized to simulate the intrinsic CMD. We produced a series of instantaneous star burst with different ages and metallicities. The ages were from the set 125 Myr, 375 Myr, 605 Myr, 855 Myr, 1.2 Gyr, 1.7 Gyr, 2.4 Gyr, 3.4 Gyr, 4.8 Gyr, 6.8 Gyr, 9.6 Gyr, 13.5 Gyr. The metallicity were from the set Z=0.003, 0.008, 0.019, 0.03, 0.05. These choices made sure that our grids were large enough to cover big age range and were small enough to discover the different features in the CMD.

The molecular clouds in the M31 bulge could complicate our investigation; we did not know the relative radial locations of individual stars, respective to the molecular clouds, so that we did not know how much foreground extinction which these stars suffered from. Fortunately, the M31 bulge is lack of molecular clouds (Li, Wang & Wakker, 2009; Melchior & Combes, 2011; Dong et al., 2016). In order to test the effects of the molecular clouds on our CMD fitting results, we performed the CMD fitting on the catalogs with or without the regions of high dust surface density. Specifically, we used the Herschel dust map (Groves et al., 2012) to identify the regions with large dust mass, i.e. the regions with 0.1 M pc. By using the dust mass to visual extinction ratio from Draine et al. (2013) and the relative extinction from Dong et al. (2014), this corresponded to A0.23 and A-A0.23. Only 24% of our FoV satisfied this requirement.

For the other parameters: 1) We adopted the distance modulus, 24.46 (780 kpc, McConnachie et al. 2005; de Grijs & Bono 2014); 2) We chose the Kroupa (2001) IMF. The exact choice of the IMF was not highly relevant because the different standard IMFs suggested in the literature differ mainly at the low stellar mass range, which lied far below our detection limit. As a result, alternative IMF could only alter the total mass of the system, but not the relative SFH; 3) The foreground Galactic extinction (A=0.282 mag, A=0.226 mag, A=0.203 mag and A=0.095 mag) was considered in the CMD fitting, while no extra local extinction was included.

We performed the CMD fitting for the WFC observation at the nine annuli defined in §3.3 and the HRC observation. We divided our fitting processes into several steps to see the contribution of stellar populations with different ages in the observed CMD. Considering that the M31 bulge is known to be dominated by old stars, we first fit the observed CMD with models older than 5 Gyr. Then, we added the stellar populations older than 1 Gyr, but younger than 5 Gyr. At the end, we used all the stellar populations.

The uncertainties of the CMD fitting include two parts: statistic and systematic errors. 1) The statistic uncertainty comes from the Poisson fluctuation of number of stars in each bin of the CMD. To determine this uncertainty, we add random errors into individual bins of the CMDs of the nine annuli and fit the new CMDs. For each annuli, this procedure has been repeated for 1000 times. The standard deviation of the output parameters are used to be the statistic uncertainty of the CMD fitting; 2) The systematic uncertainty is introduced by the distance, foreground extinction and the bin effect. The uncertainty of the distance is 0.1 mag (-0.1, 0 and 0.1 mag, de Grijs & Bono, 2014). For the local extinction, we add an extra extinction of =0, 0.06 and 0.12 into the fitting and assuming the extinction curve given by Dong et al. (2014). For the bin effect, we shift in color and magnitude with one third of the bin size (-1/3, 0 and 1/3 of the bin size). Then, we rerun the CMD fitting for these, 3 (distance) 5 (local extinction) 9 (bin shift) = 135, new dataset and derived the SFR at different ages and metallicities. The 68% percentiles of the values of these 135 new dataset are used to be the systematic uncertainty of the SFR derived above.

Figs. 15 and 16 show the final fitting results of the three steps for the WFC observation at the nine annuli with or without the high extinction regions. Overall, these two figures are very similar; That is because the high extinction regions only occupy a small fraction of our FoV. In the second column for the fitting with stellar populations older than 5 Gyr, significant residuals exist at the regions with -=2 and 23. They disappear with the inclusion of younger stellar populations. Compared the third and fourth columns, we see the potential existence of young stars with age younger than 1 Gyr, especially for the # 1 and #2 annuli.

Figs. 171819 and 20 show the results of the CMD fitting of the nine annuli: SFR as a function of look-back time; the metallicity as a function of look-back time; the cumulative mass-weighted age distribution and the mass as a function of metallicity. In all the annuli, there is a star formation activitiy around 1 Gyr ago, although it slightly shifts to the older age: This evolution is clearly shown in Fig. 21, the mass-weighted age changes from 1 Gyr in the innermost annulus to 1.5 Gyr at the outermost annulus. The relative contribution from the old stellar population increases toward the outer regions, which is due to the detection limit. In the inner region, we can see the recent star formation activity ( 500 Myr) happened in the two innermost annuli, but are not significantly beyond that. Overall, the metallicity in the M31 bulge is high for all the stars.

The top panel of Fig. 22 shows the result of the same CMD fitting for the HRC observation. However, we see that even in the last column using models with all the ages and metallicities, the model can not explain the stars with F330W-F435W1 and F435W23.5. In order to check whether these stars are blending stars, we divide the catalog into two regions using =3″ (see Fig. 2). The middle and low panels of Fig 22 shows their CMD fitting results, respectively. However, we can see that there is still many bright and red stars in the CMD outside the 3″, which can not be explained by the blending. Fig. 23 gives the produced stellar mass and the metallicity as a function of look-back time during the last 800 Myr for the whole HRC observation (black solid lines), the inner 3″ (blue dotted lines) and beyond (red dashed lines). Because the shallowness of the HRC observation, we do not trust the SFR for the stellar population older than 1 Gyr; for example, in Fig. 10, although for the extremely metal-poor metallicity (1/200 Z), we can detect the RGB stars for the 10 Gyr stellar population (red solid line), we cannot detect any stars from the 1 Gyr old stellar population (cyan dotted line). Fig. 23 tells us in the central 14″ (the FoV of HRC camera), the central 3″, and between 3″ and 14″, the star formation activity produced 1.0/0.4/0.6 and 1.8/1.0/0.6 M at around 0.375 and 0.6 Gyr ago. The stars that formed at 0.375 Gyr ago have around or sub-solar metallicity ([Fe/H]0).

4.3 Blending?

We have analyzed a very crowded field: the M31 bulge. As a result, our source catalogs could suffer from the blending effect. Compared to the HRC observation, the WFC observation is more affected by the crowding, because of the worse resolution and larger surface brightness. Furthermore, the blending is more significant at the F814W band than the F475W band, which causes stars artificially brighter and slightly redder (see also §3.4.2). Therefore the blending could potentially explain the differential CMDs between the inner six annuli and the #7 annulus shown in Fig. 14, instead of the intrinsic evolution of stellar populations along .

We used the Padova stellar isochrones and the PSF determined in §3.2 to simulate WFC observation to examine the blending effect. We scaled the SFH of the #7 annulus determined in §4.2 to let the surface densities of artificial observation match those at from 30″ to 180″ with a step size of 30″. Then, for each , we produced an artificial source catalog, multiplied the intensities of individual sources with the PSF and added them into a blank image. Poisson noise has been added into this image. After that, from this artificial image, we produced the deconvolution catalog and performed the CMD fitting. Then, we compared the output SFHs with the input one at different radii, which are shown in Fig. 24.

From this plot, we can see that even in the artificial observation with the highest surface densities, i.e. =30 arc second, the output CMD fitting does not overestimate the SFH at 2 Gyr: especially we do not see any significant star formation activities around several hundred Myr, compared to the input CMD. Therefore, the blending effect can not explain the 1 Gyr old stellar population in the central four annuli mentioned in §4.2. On the other hand, the output CMD fitting has higher SFR around 5-6 Gyr. The stars at this age range are near the detection limit. Therefore, the blending can cause dim stars with older ages ( 5 Gyr) brighter and let the CMD fitting obtain higher SFR at the 5 Gyr year old range.

5 Discussion

We first discuss whether 200 Myr old stars exist in the M31 bulge or not in §5.1. Then, we analyze the stellar populations and their spatial evolutions in the M31 bulge in §5.2. In §5.3, we present the star formation history of the M31 bulge derived from our data. At the end, we compare our results with previous ones derived from the SED fitting of photometric observation in §5.4 to examine the reliability of SED fitting on the study of stellar populations in remote galaxies.

5.1 200 Myr old stars in the M31 bulge?

Fig. 10 shows that 1) there is a handful of stars that fall into the BP region of the CMD derived from the WFC catalog; 2) a number of stars with F330W-F435W-0.3 in the HRC catalog. These stars could potentially be MS stars. However, on the other hand, evolved low-mass stars in the P-AGB, PE-AGB, AGB-manque or even the Planetary Nebulae (PNe) stage, could also be very blue, because they have lost the majority or all of their envelopes in the form of strong stellar winds. Therefore, we perform two steps below to further understand the origin of these blue stars.

First, like Rosenfield et al. (2012), we analyzed the UV CMD for the WFC stars in the BP and the stars with F435W23.3 mag in the HRC observation. The HST WFC3/UVIS F275W (2750Å) and F336W (3375Å) catalog of the PHAT survey is from Williams et al. (2014). We used a 0.1″ radius to identify the counterparts between the catalogs of Williams et al. (2014) and our ones. Fig. 25 shows the UV CMD for the stars in the BP of Fig. 10 at three different ranges defined in Fig. 12 and Fig. 26 gives the one in the HRC observation. For the WFC observation, we know that the majority of the stars with 230″, are the young massive ones within the 5 kpc star formation ring. Therefore, it is not surprising that in the CMD, these stars follow a diagonal, which can be well explained by MS stars with mass from 4 to 9 . Therefore, their age should be younger than 200 Myr. Instead, for the stars within 230″, they have a wide F275W-F336W color range. Because the extinction in the M31 bulge is very small (0.15 mag, Dong et al. 2016, which corresponds to -0.13 mag, Dong et al. 2014), the red stars can not be explained by MS stars with high extinction. Instead, these stars fall into the tracks of stars with 0.54 to 0.6 M. Therefore, these stars should be evolved low-mass stars. In the HRC observation, the stars with F330W-F435W-0.3 (blue symbols), also fall into a narrow F275W-F336W color range, like the BP stars in the WFC observation with 230″. On the other hand, the stars with F330W-F435W-0.3 have a wide color range (red symbols), like the BP stars in the WFC observation with 230″ and should be evolved low mass stars in the P-AGB, PE-AGB or AGB-manque phase.

Second, in order to further exclude the possibility that the blue stars (F330W-F435W-0.3) in the HRC observation are PNe, we cross-correlate our catalog with previous PNe catalogs. We used 1) the 12 PNe candidates in the central 21″ (80 pc) reported by Pastorello et al. (2013) using the SAURON integral-field spectroscopic observation; 2) the PNe candidates from Li Anqi, et al. (2018, in preparation). They analyzed the HST/UVIS observation from HST program, GO-12174 (Dong et al., 2014), which includes one HST/UVIS pointing (2′2′, i.e.450 pc450 pc) centered on M31* with a narrow-band filter F502N (5013Å) and a medium-band filter F547M (5475Å). 249 PNe candidates have been identified according to their unusually blue F502N-F547M colors, potentially due to the strong [O iii] 5007Å emission line, which is widely identified in PNe. We found that none of the blue stars in the HRC observation have counterparts in these two catalogs at a 0.5″ distance. Therefore, our results suggest that the blue stars in the HRC observation should not be PNe candidates, but instead should be young MS stars.

5.2 The stellar populations of the M31 bulge

In this subsection, we summarize the stellar populations in the M31 bulge.

5.2.1 Stars older than 5 Gyr

As previously known (Olsen et al., 2006; Saglia et al., 2010), the M31 bulge is dominated by stars older than 5 Gyr, which contribute more than 70% of the total mass in the innermost annulus (see Fig. 19. This is the lower limit, because of the detection limit). The old stars have super solar metallicity, with [Fe/H]0.3 at all nine annuli.

5.2.2 Stars younger than 500 Myr

Evidence of the existence of a 500 Myr old young stellar population has been shown in quantitative and qualitative studies above. First, the stars with F330W-F435W-0.3 and F435W23.3 mag have a very small F275W-F336W color range and can be explained by the tracks of young massive stars. Second, in §4.1, we mentioned that the CMD evolves along . In the inner region, there are more stars at 21F814W22 and 1F475W-F814W3, compared to the #7 annulus. These stars are not due to the blending of dim stars in the M31 bulge, as demonstrated in §4.3. Their locations in the CMD are consistent with those predicted by a stellar population with around solar metallicity and 200 Myr old. Furthermore, the SFH derived from the CMD fitting (Fig. 17 and Fig. 23) also shows that there is a star formation activity less than 500 Myr ago with roughly solar metallicity, especially the HRC observation with fine sampling and less affect from the blending. Therefore, these stars could be the relatives of young blue stars identified in the central 0.25″ by Lauer et al. (2012).

However, this stellar population only contribute 0.5% of the total stellar mass in our FoV. Fig.27 gives the mass and the mass fraction of this young stellar population. The mass fraction of young stellar population is 1.2%, then decreases to 0% in the # 5 annuli and beyond.

5.2.3 Intermediate-age Stellar population (1 Gyr)

Fig. 17 shows that in all the nine annuli, there is a SFR peak around 1 Gyr ago. The total mass of this stellar population (0.9 to 2 Gyr) is 5.4 M, 5.4% of the total stellar mass in our FoV. This stellar population contributes 9% of the total stellar mass in the innermost annuli, but decreases to 4% at # 3 annuli. The age of this stellar population evolves along too; the mass-weighted age increases from 1 Gyr to 1.5 Gyr. On the other hand, the mass-weighted metallicity decreases along . The 1 Gyr old stars can not have been contributed by the M31 disk, because no population with a similar age has been reported in previous deep HST/ACS observations in the M31 disk (Brown et al., 2006).

5.3 The SFH in the M31 bulge

Our studies above show that the star formation history in the M31 bulge is more complicated than we knew before. First, because the existence of a large number of old and metal-rich ([Fe/H]0) stars, the majority of the M31 bulge could have formed in the major merger of several metal rich systems. On the other hand, recently, Blaña Díaz et al. (2017) suggested that 2/3 of the M31 bulge could come from a box/peanut bulge due to the secular process, from their N-body simulations. Second, the star formation rate increased quickly around 1 Gyr ago in all the annuli, which suggest that this star formation activity permeated the entire M31 bulge. Star formation activities with similar ages have not been seen in the M31 disk (Brown et al., 2006). Considering the age and metallicity evolution along , this population could have arisen from secular processes; when the molecular clouds spiraled into the nuclear region, the star formation activities first happened in the outer region of the M31 bulge. The supernova explosions of massive stars contaminated the remaining molecular clouds, which kept falling into the nuclear region to trigger more star formation. Third, the potential star formation activity less than 500 Myr ago identified in the M31 nucleus is only found in the central three annuli (130″, 500 pc, see also Lauer et al. 2012), which means that this star formation event was rather localized. Instead, if the interaction between M31 and the satellite galaxy, M32, as suggested by Block et al. (2006), was real and triggered this star formation, we expect that the 200 Myr old stellar population should be distributed across a larger region which does not agree with the observation. Another possibility is that M31 used to have a central molecular zone like the one in the Milky Way. An intense star formation happened less than 500 Myr ago. Then the star formation activity depleted the molecular clouds and the feedback from the massive stars blew away the remaining gas and stopped any further star formation activity.

5.4 Comparison with the SED fitting

The M31 bulge provides us a good lab to examine the reliability to use the SED fitting of integrated stellar light to get the SFH of remote galaxies, using the output of the CMD fitting as a candle. For example, García-Benito & Pérez-Montero (2012) and Ruiz-Lara et al. (2015) test the SFH derived by both methods and they conclude that the SED fitting is a reasonable tool to derive the general properties of unresolved stellar populations.

Dong et al. (2015) used the Starburst99 to perform the SED fitting of the photometric observation from mid-UV to NIR observation of the M31 bulge. Limited by ten data points, they used only two stellar populations: one old metal-rich and another intermediate-age one. The age of the intermediate-age stellar population is 300 Myr in the central 5″. Then, it increases to 600 Myr beyond that and slowly increases to 800 Myr at 180″. The metallicity of this population is 2.5Z at the central 5″and then decrease to . The mass fraction of this stellar population is 0.2-2 percent.

Compared to Dong et al. (2015), the advantage of the current study is that we can get finer resolution in the age distribution. For example, below 2 Gyr, we find that there are two significant star formation activities: one at around 1 Gyr and the other one at less than 500 Myr ago. The luminosity-weighted age of the stellar populations 2 Gyr in our work is given as red lines in Fig. 21. In the innermost annuli (# 1), our luminosity-weighted age is 900 Myr, which is larger than 300 Myr in the central 5″, but close to the age value beyond that given in Dong et al. (2015). The difference in the central 5″ is because 1) limit by the number statistic and the detection limit of the WFC observation, the CMD fitting method cannot divide the star catalogs into small groups with the spatial resolution as high as that of the SED fitting based on the integrated light (5″); 2) the mass fraction of the 500 Myr old stellar population is low, which cannot reduce the luminosity-weighted age. However, as mentioned in §5.1, in the HRC observation, we indeed find the candidates of massive stars with F330W-F435W-0.3. The total mass of this stellar population (4 at the 375 Myr age bin) in the central 3″ derived from the CMD fitting of the HRC observation is close to 1 (in the 0-5″) given in Dong et al. (2015). The luminosity-weighted metallicity is higher than that of the SED fitting, however, which is still consistent within the large uncertainty, because in the SED fitting, we do not have any narrow-band filters which are sensitive to the metallicity variation. Therefore, in one word, the SED fitting with the integrated light can provide similar results of the CMD fitting, but with higher angular resolution and poor age and metallicity resolution. Future SED fitting of optical spectra derived from the CAHA PPAK observation of the M31 bulge can provide us with more finer age and metallicity resolution to be compared with the CMD fitting (García-Benito et al. in preparation).

6 Summary

We have investigated the stellar populations in the central 5.5′ of the M31 bulge by using HST/ACS WFC and HRC observations. Our field-of-view covers the whole M31 bulge and includes a small portion of the 5 kpc star formation ring in the north-west corner of our observation. We have produced Nyquist-sampled images and constructed source catalogs. The whole source catalog has been divided into nine annuli, centered on the M31*. We have performed the CMD fitting to study the stellar populations in the M31 bulge.

The main results have been summarized below.

  • We have demonstrated that our deconvolution method on the Nyquist-sampled image performed better in the crowded field than the pipeline described in Dalcanton et al. (2012) with the‘DOLPHOT’ software; in the artificial star cluster test, the output catalogs derived by both methods are brighter and reddened than the input stellar populations; However, the CMD produced by the deconvolution method is closer to the input one.

  • More than 2000 stars fall in the blue plume region in the WFC observation. Among them, the stars in the 5 kpc star formation ring (galactocentric radius 230″,900 pc) fall onto a diagonal line in the HST/WFC3 UVIS F275W-F336W VS F336W CMD, which indicates that they are definitely main-sequence stars. However, the stars inside the central 230″ are widely distributed in the F275W-F336W VS F336W CMD and could be explained by evolved low-mass stars, such as PAGB, PE-AGB or AGB-manque.

  • In the HRC observation of the central 14″, the stars with F435W23.3 mag could be divided into two parts by using F330W-F435W=-0.3. The F275W-F336W color of the blue one can be explained by main-sequence stars, while the red one could also be evolved low-mass stars. The former one confirms previous evidence that there is indeed an 200 Myr old intermediate-age stellar population in the M31 nucleus.

  • From the CMD fitting of the WFC and HRC observations, we can obtain the following information about the star formation history in the M31 bulge: 1) Overall, the majority of the stars in the M31 bulge, at least 70%, are old ( 5 Gyr) and metal-rich ([Fe/H]0) stars. Subsequently, the star formation rate decreased till 2 Gyr. 2) A star formation burst occurred about 1 Gyr ago throughout the entire bulge. 3) The star formation activity less than 500 Myr ago with solar metallicity was only concentrated in the central 130″ (500 pc).

  • The identification of the 1 Gyr old star formation activity confirms our previous results derived from the spectral energy distribution fitting of the photometric observation at ten bands of the HST WFC3/ACS observation. The radial gradients of age and metallicity of these stars suggest that they may have formed during the infall of molecular clouds from the M31 disk (Opitsch et al., 2018).

  • The spatial distribution of the stellar population less than 500 Myr old is inconsistent with the scenario that the potential collision between the M31 and M32 triggered the star formation, which should produce star formation in a larger scale. Instead, we suggest that they could be due to the nuclear star formation, like the central molecular zone in the Milky Way. Then the molecular cloud has been consumed by star formation activity and also been destroyed by the feedback.

  • The properties of the stellar populations in the M31 bulge derived by the spectral energy distribution fitting of the integrated light is similar to those from the CMD fitting, which suggests that spectral energy distribution fitting method is still a reasonable method to study the stellar populations in remote galaxies.


The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n [614922]. and it was also supported by NASA via the grant GO-12055, provided by the Space Telescope Science Institute. This work uses observations made with the NASA/ESA Hubble Space Telescope and the data archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. H. D. acknowledges the hospitality of University of Massachusetts, Amherst during his visit. We are grateful to Q. Daniel, Wang and Philip Rosenfield for many valuable comments and discussion and Pauline Barmby for providing the Spitzer/IRAC image.


  • Barmby et al. (2006) Barmby, P., Ashby, M. L. N., Bianchi, L., et al. 2006, ApJ, 650, L45
  • Beaton et al. (2007) Beaton, R. L., Majewski, S. R., Guhathakurta, P., et al. 2007, ApJ, 658, L91
  • Bender et al. (2005) Bender, R., Kormendy, J., Bower, G., et al. 2005, ApJ, 631, 280
  • Blaña Díaz et al. (2017) Blaña Díaz, M., Wegg, C., Gerhard, O., et al. 2017, MNRAS, 466, 4279
  • Block et al. (2006) Block, D. L., Bournaud, F., Combes, F., et al. 2006, Nature, 443, 832
  • Brown et al. (1998) Brown, T. M., Ferguson, H. C., Stanford, S. A., & Deharveng, J.-M. 1998, ApJ, 504, 113
  • Brown et al. (2006) Brown, T. M., Smith, E., Ferguson, H. C., et al. 2006, ApJ, 652, 323
  • Chabrier (2003) Chabrier, G. 2003, ApJ, 586, L133
  • Courteau et al. (2011) Courteau, S., Widrow, L. M., McDonald, M., et al. 2011, ApJ, 739, 20
  • Conroy & van Dokkum (2016) Conroy, C., & van Dokkum, P. G. 2016, ApJ, 827, 9
  • Crane et al. (1992) Crane, P. C., Dickel, J. R., & Cowan, J. J. 1992, ApJ, 390, L9
  • Dalcanton et al. (2009) Dalcanton, J. J., Williams, B. F., Seth, A. C., et al. 2009, ApJS, 183, 67
  • Dalcanton et al. (2012) Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18
  • Davidge (2001) Davidge, T. J. 2001, AJ, 122, 1386
  • Davidge et al. (2005) Davidge, T. J., Olsen, K. A. G., Blum, R., Stephens, A. W., & Rigaut, F. 2005, AJ, 129, 201
  • Depoy et al. (1993) Depoy, D. L., Terndrup, D. M., Frogel, J. A., Atwood, B., & Blum, R. 1993, AJ, 105, 2121
  • Dong et al. (2011) Dong, H., Wang, Q. D., Cotera, A., et al. 2011, MNRAS, 417, 114
  • Dong et al. (2014) Dong, H., Li, Z., Wang, Q. D., et al. 2014, ApJ, 785, 136
  • Dong et al. (2015) Dong, H., Li, Z., Wang, Q. D., et al. 2015, MNRAS, 451, 4126
  • Dong et al. (2016) Dong, H., Li, Z., Wang, Q. D., et al. 2016, MNRAS, 459, 2262
  • Draine et al. (2013) Draine, B. T., Aniano, G., Krause, O., et al. 2013, arXiv:1306.2304
  • Dressler & Richstone (1988) Dressler, A., & Richstone, D. O. 1988, ApJ, 324, 701
  • Fisher & Drory (2010) Fisher, D. B., & Drory, N. 2010, ApJ, 716, 942
  • Garcia et al. (2010) Garcia, M. R., Hextall, R., Baganoff, F. K., et al. 2010, ApJ, 710, 755
  • García-Benito & Pérez-Montero (2012) García-Benito, R., & Pérez-Montero, E. 2012, MNRAS, 423, 406
  • Groves et al. (2012) Groves, B., Krause, O., Sandstrom, K., et al. 2012, MNRAS, 426, 892
  • Girardi et al. (2012) Girardi, L., Barbieri, M., Groenewegen, M. A. T., et al. 2012, Astrophysics and Space Science Proceedings, 26, 165
  • de Grijs & Bono (2014) de Grijs, R., & Bono, G. 2014, AJ, 148, 17
  • King et al. (1995) King, I. R., Stanford, S. A., & Crane, P. 1995, AJ, 109, 164
  • Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  • Kormendy (1988) Kormendy, J. 1988, ApJ, 325, 128
  • Kormendy & Kennicutt (2004) Kormendy, J., & Kennicutt, R. C., Jr. 2004, ARA&A, 42, 603
  • Lauer et al. (1998) Lauer, T. R., Faber, S. M., Ajhar, E. A., Grillmair, C. J., & Scowen, P. A. 1998, AJ, 116, 2263
  • Lauer (1999) Lauer, T. R. 1999, PASP, 111, 227
  • Lauer et al. (2012) Lauer, T. R., Bender, R., Kormendy, J., Rosenfield, P., & Green, R. F. 2012, ApJ, 745, 121
  • Lewis et al. (2015) Lewis, A. R., Dolphin, A. E., Dalcanton, J. J., et al. 2015, ApJ, 805, 183
  • Li, Wang & Wakker (2009) Li, Z., Wang, Q. D., & Wakker, B. P. 2009, MNRAS, 397, 148
  • Li et al. (2011) Li, Z., Garcia, M. R., Forman, W. R., et al. 2011, ApJ, 728, L10
  • Marigo et al. (2013) Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488
  • Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77
  • McConnachie et al. (2005) McConnachie, A. W., Irwin, M. J., Ferguson, A. M. N., et al. 2005, MNRAS, 356, 979
  • McGaugh & Schombert (2014) McGaugh, S. S., & Schombert, J. M. 2014, AJ, 148, 77
  • Meidt et al. (2014) Meidt, S. E., Schinnerer, E., van de Ven, G., et al. 2014, ApJ, 788, 144
  • Melchior & Combes (2011) Melchior, A.-L., & Combes, F. 2011, A&A, 536, A52
  • Melchior & Combes (2013) Melchior, A.-L., & Combes, F. 2013, A&A, 549, A27
  • Monachesi et al. (2011) Monachesi, A., Trager, S. C., Lauer, T. R., et al. 2011, ApJ, 727, 55
  • Olsen et al. (2006) Olsen, K. A. G., Blum, R. D., Stephens, A. W., et al. 2006, AJ, 132, 271
  • Opitsch et al. (2018) Opitsch, M., Fabricius, M. H., Saglia, R. P., et al. 2018, A&A, 611, A38
  • Pastorello et al. (2013) Pastorello, N., Sarzi, M., Cappellari, M., et al. 2013, MNRAS, 430, 1219
  • Rich et al. (1993) Rich, R. M., Mould, J. R., & Graham, J. R. 1993, AJ, 106, 2252
  • Rich & Mighell (1995) Rich, R. M., & Mighell, K. J. 1995, ApJ, 439, 145
  • Rosenfield et al. (2012) Rosenfield, P., Johnson, L. C., Girardi, L., et al. 2012, ApJ, 755, 131
  • Rosenfield et al. (2016) Rosenfield, P., Marigo, P., Girardi, L., et al. 2016, ApJ, 822, 73
  • Ruiz-Lara et al. (2015) Ruiz-Lara, T., Pérez, I., Gallart, C., et al. 2015, A&A, 583, A60
  • Saglia et al. (2010) Saglia, R. P., Fabricius, M., Bender, R., et al. 2010, A&A, 509, A61
  • Sandage (1990) Sandage, A. 1990, JRASC, 84, 70
  • Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
  • Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., Weinberg, M. D., Schneider, S., Carpenter, J. M., Beichman, C., Capps, R., Chester, T., Elias, J. et al., 2006, AJ, 131, 1163S
  • Stephens et al. (2003) Stephens, A. W., Frogel, J. A., DePoy, D. L., et al. 2003, AJ, 125, 2473
  • Toomre (1977) Toomre, A. 1977, Evolution of Galaxies and Stellar Populations, 401
  • Vázquez & Leitherer (2005) Vázquez, G. A., & Leitherer, C. 2005, ApJ, 621, 695
  • Williams et al. (2014) Williams, B. F., Lang, D., Dalcanton, J. J., et al. 2014, ApJS, 215, 9

Figure 1: (a): GALEX composite mosaic of M31 (red: NUV, blue: FUV). The polygons show the boundary of Brick 1; Green and cyan polygons are halves of Brick 1 with and without half-pixel dithering strategy, respectively. (b): HST ACS/WFC F814W mosaic. The units of the mosaic are counts . The green boxes and labels represent individual ‘Fields’. The cyan polygon shows the field-of-view of the ACS/HRC observation.

Figure 2: The ACS/HRC F330W image of the M31 nucleus (29″26″). The units of the image are counts . The radius of the green circle is 3″.

Figure 3: The recovered fraction as a function of input magnitude derived from the artificial star tests for the WFC (top panels) and HRC (bottom panels) observations. In the top panel, the solid lines from left to right represent the annuli from innermost to outermost with decreasing surface densities.

Figure 4: The difference between the output and input magnitudes at the F475W (top row) and F814W (middle row) bands for the innermost and outermost annuli of the WFC observation, as well as the HRC observation (bottom row). The red lines and bars represent the mean and standard deviation of the differences at individual input magnitude bin. The red dotted and dashed lines represent the 70% and 50% completeness limits.

Figure 5: The Hess diagrams for Field 11 (top panel) and Field 18 (bottom panel), with the DOLPHOT ‘st’ catalogs from PHAT pipeline team (left column) and our deconvolution catalogs (right column). These two ‘fields’ have the highest and lowest star number surface densities in our study. Compared to Field 11, the detection limit of Field 18 has increased by 1 mag at the F814W band. The black solid lines in the Hess diagram of these two ‘Fields’ roughly mark the detection limit of the DOLPHOT ‘st’ catalog. Our deconvolution method finds more dim and red stars in the low right part of the Hess diagram. The sharp cut of the Hess diagram of our deconvolution catalogs at the dim and red side is due to the S/N cut at both bands. Compared to our deconvolution catalog, the DOLPHOT ‘st’ catalog has extra sources at 1) F475W-F814W6 and 2)F475W-F814W0 and F814W27, which are due to their S/N cut at only one band.

Figure 6: The LF for Field 11 (top panel) and Field 18 (bottom panel), at F475W (left column) and F814W (right column) with the DOLPHOT ‘st’ catalogs from PHAT pipeline team (dashed lines) and our deconvolution catalogs (solid lines). Our method detects more dimmer sources. The shape break of the LF from our deconvolution catalog is due to the S/N cut at both bands. Instead, the long tail in the LF from the DOLPHOT ‘st’ catalog is also due to their S/N cut at only one band.

Figure 7: Comparison of the CMD of DOLPHOT (Middle column) and deconvolved catalog (Right column) in four selected artificial star clusters. Left column is the WFC F814W images of individual star clusters. The radii of the red and magenta circles are the effective radii and 2”. We extract the CMD in the middle and right columns from the magenta circles. In the middle and right columns, the black circles are the input artificial sources, while the blue and magenta circles are the sources detected by DOLPHOT and deconvolution methods, respectively.

Figure 8: The difference between the input and output parameters; magnitude and color, as a function of the input magnitude for two artificial star clusters.The red lines represent the same input and output parameters. The three cyan lines are the mean and one standard deviation of the difference between the input and output parameter. In some panels, the cyan lines overlays on the red lines. The data in the first and third rows are from the DOLPHOT catalog, while those in the second and fourth rows are from the deconvolution catalog.

Figure 9: The surface brightness (mag arcsec) at the F814W band along the galactocentric radius (, in units of arc second). The region inside the 230″ has brighter surface brightness than that of artificial cluster ID 1, which indicates that the deconvolution catalog should be preferred in our study.

Figure 10: The Hess diagram of the M31 bulge from the WFC (top) and HRC (bottom) observations. On the top left panel, the cyan line marks the saturation limit of the WFC observation, while the blue and red parallelograms embed blue and red stars, which are discussed in §4.1. The green lines in the left column are the 50% detection limit. We overlay the isochrones from the Padova stellar evolutionary tracks in the right column. The solid and dashed lines represent the stellar populations with 1/200 and 2.5 solar metallicities, respectively, at five different ages, 0.1 (purple), 0.5 (blue), 1 (cyan), 5 (green), 10 Gyr (red).

Figure 11: The Hess diagram of the stars in the M31 bulge from the WFC observation at the nine annuli defined in Table 2. Only the stars brighter than 50% completeness limit at both bands are kept. We also overlay the isochrones from the Padova stellar evolutionary tracks. The solid and dashed lines represent the stellar populations with 1/200 and 2.5 solar metallicities, respectively, at five different ages, 0.1 (purple), 0.5 (blue), 1 (cyan), 5 (green), 10 Gyr (red).

Figure 12: The spatial distribution (left) and surface number density (right) of sources within the blue parallelogram in the top left panel of Fig. 10. The two circles have galactocentric radii of 80″ and 230″, respectively. For each galactocentric radius, the surface number density is the ratio of the number of stars in an annulus with 5″ width to the effective sky area (arc second) within this annulus.

Figure 13: The CMD of stars within the blue parallelogram in the top left panel of Fig. 10, but at different galactocentric radii.

Figure 14: The difference of the Hess diagrams between the six inner annuli and the #7. The recovered fraction between the annuli is corrected. The CMD of the #7 annulus have been scaled to have the same total number of stars as the other annuli. The green and the red lines are the 200 Myr and 500 Myr old stellar populations with solar metallicity. The cyan box encompasses a region in the CMD of the #1 annulus with significantly extra stars relative to the #7 annulus.

Figure 15: The CMD fitting result for the WFC observation. The first column is the observed CMD. The second, third and fourth columns are the fitting residual using the models with 5 Gyr only, 1 Gyr only and all the ages.

Figure 16: Similar to Fig. 15, but for the regions without high extinction (see §4.2).

Figure 17: The star formation history for the nine annuli derived from the CMD fitting with (solid lines) or without high-extinction regions (red diamonds, see §4.2).

Figure 18: The metallicity evolution for the nine annuli derived from the CMD fitting with (solid lines) or without high-extinction regions (red diamonds, see §4.2).

Figure 19: The cumulative mass-weighted age distribution for the nine annuli derived from the CMD fitting with (solid lines) or without high-extinction regions (red diamonds, see §4.2).

Figure 20: The mass as a function of metallicity for the nine annuli derived from the CMD fitting with (solid lines) or without high-extinction regions (red diamonds, see §4.2).

Figure 21: The evolution of the mass-weighted (black lines) and luminosity-weighted (red lines) age (left) and metallicity (right) for the stellar populations younger than 2 Gyr at the nine annuli. The black/red solid and dashed lines are for the CMD fitting with and without the high-extinction regions (see §4.2). The error bars include both statistic and systematic errors in the CMD fitting described in §4.2. The blue solid and dashed lines are the age/metallicity and their uncertainties of the intermediate-age stellar population given in Dong et al. (2015). The eight black vertical dotted lines are the boundaries of the nine annuli defined in Table 2.

Figure 22: The CMD fitting result for the HRC observation. The first column is the observed CMD. The second, third and fourth columns are the fitting residual using the models with 5 Gyr only, 1 Gyr only and all the ages.

Figure 23: Left: total stellar mass as a function of look-up time up to 0.8 Gyr; Right: age-metallicity distribution. The black, blue and red lines are for all the stars, the stars inside and outside =3″ of the HRC observation.

Figure 24: The CMD fitting outputs from the artificial observation with the surface densities equal to those at different radii. The solid line is the input SFH, while the dotted line is the output SFH.


Figure 25: The CMD (F275W-F336W vs F336W) of stars fall in the blue parallelogram of Fig. 10 at three different regions (from top to bottom panels) defined in Fig. 12. In the left and right columns, we compared the CMD with the MS tracks for young stars with 9, 7, 5, 4 and 3 solar mass (left) and the tracks of the evolved low-mass stars at the P-AGB (green lines, from top to bottom, 0.9, 0.75 and 0.6 ) and HP-HB (red lines, from top to bottom, 0.58, 0.56, 0.54, 0.5 and 0.48 .) phases (right), as described in Rosenfield et al. (2012).

Figure 26: The CMD (F275W-F336W vs F336W) of stars detected in the HRC observation with F336W23.3 mag. The blue and red symbols represent the stars with F336W-F435W-0.3 and F336W-F435W-0.3 respectively, while the diamonds and triangles inside and outside the galactocentric radius of 3″. In the left and right columns, we compared the CMD with the MS tracks for young stars with 9, 7, 5, 4 and 3 solar mass (left) and the tracks of the evolved low-mass stars at the P-AGB (green lines, from top to bottom, 0.9, 0.75 and 0.6 ) and HP-HB (red lines, from top to bottom, 0.58, 0.56, 0.54, 0.5 and 0.48 .) phases (right), as described in Rosenfield et al. (2012).

Figure 27: The total mass (black lines) and mass fraction (red lines) of the stellar population of 500 Myr old (dotted lines) and 900 MyrAge 2 Gyr old (solid lines). The eight black vertical dotted lines are the boundaries of the nine annuli defined in Table 2.
Filter Central Wavelength (Å) Width (Å) Exposure
F475W 4760 1458 10+360+360+470+700
F814W 8333 2511 15+350+800+550
F330W 3354 588 8298
F435W 4297 1038 6724+6798

Note. – : the solid numbers mean that the corresponding dithered exposures are used to produce the Nyquist Images for the WFC observations in §3.

Table 1: Broad-band Filters
ID Radius F475W(70%) F475W(50%) F814W(70%) F814W(50%)
1 0″-69″ 24.2 24.6 21.9 22.3
2 69″-101″ 24.6 25.0 22.4 22.8
3 101″-126″ 24.8 25.2 22.7 23.1
4 126″-147″ 25.0 25.4 22.9 23.3
5 147″-169″ 25.2 25.5 23.1 23.5
6 169″-193″ 25.3 25.7 23.2 23.7
7 193″-218″ 25.4 25.8 23.4 23.9
8 218″-249″ 25.6 26.0 23.6 24.0
9 249″-339″ 25.9 26.3 23.9 24.4

Note. – The 70% and 50% completeness limit of the WFC observation at nine annuli of the F475W and F814W bands.

Table 2: Completeness limit
log Age Total Mass Radius MD SB
ID (yr) (M) (pc) (M pc) (mag arcsec) (mag)
1 9.6 9.6e+04 3.2 3055.8 18.0 0.2
2 7.8 2.7e+03 3.0 95.0 19.6 0.3
3 7.8 3.5e+03 2.7 150.6 19.0 0.6
4 9.4 7.5e+03 1.5 1126.8 20.1 1.5
5 8.1 6.2e+02 3.2 19.3 21.0 0.5
6 9.8 9.7e+03 3.1 313.4 19.6 0.3
7 6.7 8.6e+02 2.7 38.4 20.0 0.4
8 8.4 1.8e+03 8.8 7.3 20.6 1.6
9 8.1 5.4e+02 1.9 45.0 20.7 0.3
10 7.3 5.8e+02 2.7 25.5 20.6 1.8
11 9.6 1.8e+04 5.8 171.7 19.4 0.9
12 8.6 1.2e+03 2.3 74.3 20.4 0.2
13 9.7 9.3e+03 4.7 133.7 21.0 0.8
14 9.9 5.0e+03 2.7 219.0 21.0 1.3
15 9.6 2.0e+04 1.3 3926.0 19.7 2.1
16 9.6 4.0e+04 3.7 927.5 19.1 0.3
17 9.9 6.5e+04 2.8 2699.5 20.8 2.9
18 8.1 3.5e+02 3.1 11.4 20.5 1.6
19 9.7 9.1e+04 1.8 8511.0 20.9 2.8
20 7.1 1.1e+03 7.0 7.4 20.0 0.9
21 7.8 3.5e+03 1.6 442.6 20.9 1.4
22 7.4 2.5e+02 1.4 40.2 19.7 1.3
23 8.2 1.8e+03 3.6 43.4 20.6 0.8
24 8.6 1.4e+04 1.9 1198.3 18.2 1.1
25 8.6 7.7e+02 1.7 87.3 20.6 1.3
26 6.8 3.3e+02 1.8 31.6 20.0 1.4
27 8.8 2.4e+03 1.6 278.4 20.2 1.4
28 8.6 7.8e+03 2.8 307.2 20.9 0.3

Note. – : the effective radius. : the mass surface density within the effective radius. : the F814W surface brightness within the central 0.5″ (1.9 pc).

Table 3: Artificial Star Clusters
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description