A GLOBAL STAR-FORMING EPISODE IN M31 2–4 GYR AGO

Historically it has been challenging to model the stellar populations of M31 inside of 20 kpc because of the high dust content which causes significant differential reddening (e.g., Dalcanton et al. 2012; Draine et al. 2014). To overcome this difficulty, we take advantage of the large contiguous area of the PHAT survey to find small regions containing very little extinction.

The photometry sample for our study was chosen using the extinction maps of Dalcanton et al. (2015), where dozens of 15'' × 15'' (60 × 250 pc, deprojected) dust-free regions were found across the PHAT footprint by measuring the width of the red giant branch (RGB) in the NIR. We selected 365 of these regions, and grouped them into eight distinct areas to provide a wide range of radius and azimuth: five areas along the major axis, one along the minor axis, and two areas of large deprojected radius.

Figure 1 shows our groups plotted on both a map of the RGB width within the PHAT footprint as well as a 3.6 micron Spitzer map of M31. Table 1 provides the median, minimum, and maximum deprojected radii, total area, total number of stars, and best-fit differential reddening for each group. The group numbers in the table correspond to the labels on the figure.

Zoom In Zoom Out Reset image size

Table 1.  Properties of the Eight Samples Taken from the PHAT Dust-free Regions

Region	Rdeproj (kpc)	Mindeproj (kpc)	Maxdeproj (kpc)	Area (kpc2)	N$_{{\rm stars}}$	Σ (arcsec−2)	dAV (mag)
1	3.00	2.69	3.46	0.075	122561	22.7	0.6
2	4.69	4.13	4.93	0.034	50177	20.3	0.6
3	7.81	7.19	8.51	0.182	198732	15.2	0.4
4	14.24	13.42	14.92	0.075	47038	8.7	0.2
5	16.94	15.32	19.35	0.245	131251	7.5	0.2
6	17.67	17.17	18.52	0.044	15078	4.8	0.0
7	18.34	17.19	19.87	0.288	91090	4.3	0.0
8	18.81	17.15	19.93	0.191	48055	3.5	0.2

Download table as:  ASCIITypeset image

The CMDs of six of our eight groups, which span the range of photometric quality of our data, are shown in Figure 2. While the RGB of samples from smaller radii are broader due to crowding effects (see Williams et al. 2014, for details), these RGBs are narrow and well-defined. In addition, the red clump feature, visible at F814W∼24.5 in all CMDs sensitive to that depth, shows very little, if any, extension. These CMDs are consistent with the Dalcanton et al. (2015) measurement of these regions containing <0.6 mag of differential extinction. As a result, we allowed up to 0.6 mag of differential extinction (dA_V) in our model fits. Including this small amount of differential extinction improved the fits, and did not change the measured age distributions beyond our measured uncertainties. However, it did impact some of the details of the SFHs at intermediate ages in the two innermost groups. These groups were best-fit with dA_V = 0.6 mag. The exact choice of dA_V had no effect on our result within our uncertainties. However, if no differential extinction was included at these radii, the results differed significantly, and were considerably less consistent with the results at larger radii.

Zoom In Zoom Out Reset image size

We chose these regions because of their lack of dust. Therefore, they are biased against containing recent star formation. Thus, we have confined our comparison across the disk to ages >500 Myr. Ages younger than this are explored in detail in Lewis et al. (2015) and are highly structured spatially. By 500 Myr, the structure is still apparent, though less pronounced. In any case, this bias is unavoidable for our present analysis, so we will take it into consideration in any interpretations of our results.

Each of these groups has an appropriate set of artificial stars measured from a region of similar stellar density. To maximize our age sensitivity, we fit the optical CMDs of each set of photometry and artificial stars using the software package MATCH, using techniques well documented in the literature (Dolphin 2002, 2012, 2013; Gallart et al. 2005; Weisz et al. 2014). In short, the CMD is binned into a 2D histogram, and this histogram is fitted by a linear combination of model histograms from a fine grid of age and metallicity. The best-fit provides the most likely age and metallicity distribution for the stars. We then apply a Monte Carlo technique for determining the range of SFHs that would provide acceptable fits to the data. The results of the Monte Carlo runs yield random uncertainties for our measurements given the reliability of the models (Dolphin 2013).

To determine the significance of a particular feature in the best-fitting SFH, one must account for systematic uncertainties that come from deficiencies in the models or offsets between the models and measurements. The technique for determining these uncertainties is detailed in Dolphin (2012). Briefly, the photometry is fitted to the models in 50 independent runs. In each run, the models are shifted with respect to the data by a random and small amount (±0.02 in log of the effective temperature and ±0.17 in bolometric magnitude). The results of the 50 runs are merged to calculate the variance of the results in each time bin, which is adopted as the systematic uncertainty. These uncertainties were combined with the random uncertainties to produce our total uncertainties on our age distributions.

We fit our data using several different stellar evolution models. We adopted fits with the Padova models (Girardi et al. 2000; Marigo et al. 2008) as our fiducial SFHs; however, we show the nominal SFHs using the BASTI (Pietrinferni et al. 2004) and PARSEC (Bressan et al. 2012) stellar evolution libraries to give a sense of the effects that different models have on age when fitting data of our varying depths.

We also produced model CMDs of constant star formation rate over the past 14 Gyr, convolved with the errors, completeness, and reddening measured from the PHAT data of varying depths. We fitted these using the same techniques used to fit the observed CMDs. As with the extensive tests of such fitting in the literature (e.g., Dolphin 2002), in our tests, the measured SFHs fell within the uncertainties of the input constant SFH at the full range of photometric depths included in our data.

In addition to our new measurements, we tested our ability to recover reliable SFHs using data of the depth and quality of the PHAT photometry. Our test data were taken from the HST archive using the ACS field in the warp in the outer regions of the southern M31 disk, previously analyzed in detail by Bernard et al. (2012). This field has very deep data and contains a clear feature at ∼2 Gyr in the published age distribution. We analyzed this field in two ways. First, we ran the entire data set through our entire photometry and CMD-fitting technique with the Padova models to ensure that we could reproduce the Bernard et al. (2012) result. Then, we reduced a subset of the data to mimic the depth of the PHAT photometry in order to look for biases in the distribution of ages younger than 5 Gyr that could potentially be attributed to the shallower nature of the PHAT survey. Even using only a small subset of the Bernard et al. (2012) data resulted in relatively deep photometry due to the low stellar density of the far outer disk, so we limited the fitting to a portion of the CMD similar to the portion we were able to fit using the PHAT data in our inner fields (see black line in the upper-right panel of Figure 3). Thus, when fitting the shallow data, the fitting routine had no access to the subgiant branch or lower main-sequence, as these features fall outside of the included portion of the shallow CMD.

Zoom In Zoom Out Reset image size

Figure 3 shows comparison between the age distributions younger than 5 Gyr from fitting the full depth of the data and from the much shallower subset of the data. We do not include ages greater than 5 Gyr on the plots to focus on the performance in the epoch of interest.¹¹ Thus the cumulative fractions are relative to only the subset of stars with ages <5 Gyr, not the total stellar population. We outline the timing of the burst we find in the inner disk by shading different epochs. The age distribution from the full depth is shown with the dashed line. As published in Bernard et al. (2012) this SFH has a lull in star formation at 4–5 Gyr followed by a sharp burst at 2–3 Gyr. The fit to the shallower data also detects a burst at 2–3 Gyr, but ascribes ∼25% more of the <5 Gyr old stellar population to the burst than the fit to the deep data. Thus, both measurements agree on the presence of the burst, though its exact strength is more uncertain for the shallower data. However, the two measurements are fully consistent within the uncertainties, shown as the gray shaded region, confirming that our technique accurately estimates the precision with which we can determine the SFH in this age range.

Fits to all model sets show that the population of stars with ages <5 Gyr is dominated by a strong episode of star formation. All of our SFHs show that star formation essentially shuts down over the past 1–2 Gyr, after undergoing a significant burst of star-forming activity 2.0–3.5 Gyr ago. The possible exception is in the group near 7.8 kpc, where star formation was more steady.

Since our data are not of sufficient depth to probe the main sequence turnoff for stars of ages 2–3 Gyr, it is not necessarily instantly obvious which CMD feature is driving the detection of the large population of stars with these ages in our data. To qualitatively assess the detection, we produced model CMDs using the best fit and omitting the star formation from 2.0–3.5 Gyr, using the 4.7 kpc group as a fiducial. These model CMDs, along with the observed CMD and the difference between the two, are shown in Figure 5. The driving features appear to be the bright-blue portion of the red clump, the asymptotic giant branch, and the relatively blue and vertically oriented portion of the RGB. These features are apparently fit best by models of 2.0–3.5 Gyr in age. We have verified that models of older age with lower metallicity provide a poorer fit, as do models of constant star formation rate.

Zoom In Zoom Out Reset image size

We note that these CMD features are not the most well-tested in the stellar evolution models (e.g., Gallart et al. 2005). Therefore, there is a chance that future changes to models of the RGB, asymptotic giant branch, and red clump and/or deeper data could prove the population consistent with a more constant SFH. However, the detection of a burst of the same age in 14 deeper fields in the outer disk (25–90 kpc deprojected radii) at a variety of azimuthal angles (Bernard et al. 2015), along with the agreement between our fields of depths varying by two magnitudes and covering 17 kpc of the inner disk in deprojected galactocentric distance, strongly suggests that our result is robust given the models currently available for fitting resolved stellar photometry.

3.1.1. Age of the Burst

3.1.2. Stellar Mass Produced

To investigate the amount of stellar mass formed in this burst, we calculated the surface density of the stellar mass formed from 2.0 to 3.5 Gyr ago as a function of median radius, and we compared this with the surface density of stars formed over other epochs of similar duration from 0.5–5.0 Gyr ago. These epochs are marked on Figure 4. We avoid the most recent 500 Myr as our regions were chosen to avoid dust, making us biased against recent star formation, and because this time period is investigated in detail in Lewis et al. (2015). We plot these profiles in the left panel of Figure 6. We note that the results from the shallower data at small radii are consistent with an extrapolation of the measurements from deeper data at larger radii. The right panel shows a sand-pile histogram of the fraction of the total 0.5–5 Gyr old stellar mass formed in each epoch.

Zoom In Zoom Out Reset image size

Integrating the exponential profile plotted for the 2.0–3.5 Gyr time period shown in green in Figure 6 (see Table 2), assuming azimuthal symmetry, yields 8.1 × 10⁹ ${{M}_{\odot }}$ of stars formed. Integrating the exponential profile for the 3.5-5.0 Gyr period yields only 2.5 × 10⁹ ${{M}_{\odot }}$, and the profile for the 0.5–2.0 Gyr epoch also yields only 2.2 × 10⁹ ${{M}_{\odot }}$. Thus, this was a galaxy-wide factor of ∼3–4 increase in the star formation rate averaged over 1.5 Gyr, resulting in the production of 6 × 10⁹ ${{M}_{\odot }}$ of additional stellar mass. Assuming ∼30% star formation efficiency, this mass would be associated with 2 × 10¹⁰ ${{M}_{\odot }}$ of baryonic matter.

Table 2.  Exponential Profile Fit Parameters and Integrated Stellar Masses for the Three Age Ranges

Age Range	Scale Length (kpc)	Normalization (${{M}_{\odot }}$ kpc−2)	Integrated Stellar Mass (${{M}_{\odot }}$)
0.5–2.0 Gyr	7.1 ± 1.0	8.9$\;\pm \;3.7\times {{10}^{6}}$	2.2 ± 1.1 × 109
2.0–3.5 Gyr	4.5 ± 0.5	6.6$\;\pm \;2.1\times {{10}^{7}}$	8.0 ± 2.5 × 109
3.5–5.0 Gyr	4.9 ± 0.8	1.8$\;\pm \;1.0\times {{10}^{7}}$	2.5 ± 1.1 × 109

Download table as:  ASCIITypeset image

This large amount of star formation likely significantly depleted the gas supply, resulting in the low SFR seen today. The current gas mass in M31 is ∼7 × 10⁹ ${{M}_{\odot }}$ (Draine et al. 2014), and the star formation rate is ∼0.7 ${{M}_{\odot }}$ yr⁻¹, yielding a depletion timescale of ∼10 Gyr. However, the low SFR of M31 puts it into a transition zone between the blue and red bimodal color populations of large galaxy surveys (Mutch et al. 2011). Thus, this recent burst may mark an event in M31 history that triggered the current global transformation toward the red sequence. We consider whether this star formation 2.0–3.5 Gyr ago may signify an interaction or merger below.

It is possible that an interaction with another Local Group galaxy triggered the star formation episode in M31 2–4 Gyr ago. If so, based on previous studies, the most likely candidates are M33 and/or M32.

The timing of the burst at 2–3 Gyr ago is consistent with that seen in the outer M31 disk by Bernard et al. (2015) with deeper data over 14 fields spread all over the outer M31 disk. Simulations that reproduce the current velocities and positions of M31 and M33 suggest a close passage 2–3 Gyr ago(McConnachie et al. 2009). Furthermore, a burst of similar age is seen in the star formation history of M33 (see SFHs in Williams et al. 2009). Numerical simulations do show that interactions can increase star formation rates by factors of four to five (e.g., Springel 2000), but such increases are only seen in simulations of very close interactions between galaxies of similar mass.

There are some characteristics of the M31–M33 system that make the M31–M33 interaction explanation for such a star formation episode less clear. First is the amount of mass in stars formed. Second is galaxy mass and morphology.

Our measurements suggest that this burst was widespread and produced a large amount of stellar mass in the massive galaxy M31. A fly-by of a galaxy ∼10% of the mass of M31 may not be provide enough of a perturbation and may not include enough gas to explain this amount of increased star formation.

Furthermore, the masses and morphologies of the two galaxies do not seem consistent with a strong interaction between them. When smaller galaxies interact with larger galaxies, typically the morphology of the less massive galaxy is more affected than that of the more massive one (e.g., Mayer et al. 2001; Governato et al. 2004). M33 has a relatively undisturbed morphology, showing only a warped outer disk and a disk break beyond 8 kpc (Ferguson et al. 2007; McConnachie et al. 2009; Williams et al. 2009) despite being significantly less massive. M31 has a highly evolved morphology, as if it has been subjected to strong interactions, forming a large bulge and relatively quiescent disk. Therefore, the morphologies of the two galaxies (especially that of M33) are not those expected if they have recently been through a strong interaction.

Still, the simultaneous nature of the episodes seen in M31 and M33 suggest that there could have been an interaction, and the simulations of McConnachie et al. (2009) suggest that such an interaction occurred ∼2–3 Gyr ago. Perhaps future simulations will show that a major star formation episode can be triggered in a massive galaxy without causing major morphological changes to a low-mass galaxy. Perhaps such a burst could occur if the M31 disk was relatively gas-rich at the time and if the interaction had just the right impact parameter and orientation to induce some kind of resonance in the disk. Detailed hydrodynamic simulations of all possible interactions between the systems may be required to answer this question.

On the other hand, the M31–M32 system appears more consistent with having undergone a strong interaction a few Gyr ago. Block et al. (2006) simulate a collision between M31 and M32 ∼200 Myr ago, which is more recent than our measurements probe, but suggests that the system has interacted strongly in the past. Morphologically, as a very compact dwarf elliptical very close to M31, M32 appears more likely to have been strongly influenced by an interaction with M31 in the recent past. Furthermore, there is a significant 2–5 Gyr old population in M32 which comprises 40% of its mass (Davidge & Jensen 2007; Monachesi et al. 2012). It is possible that M32 is the remnant of a galaxy that was once more massive and gas rich (Bekki et al. 2001). This relatively gas rich galaxy could have strongly interacted with M31 about 2 or 4 Gyr ago, triggering the burst in M31 as well as the centrally concentrated star formation in M32, generating the 2–5 Gyr old population in M32, and leaving M32 the compact dwarf elliptical galaxy we see today. On the other hand, M32 does not currently show kinematic or structural evidence of a strong tidal disturbance (Lauer et al. 1998; Howley et al. 2013).

There is also the possibility that both of these scenarios occurred during the epoch in question. Perhaps M33 was much closer to M31 during this epoch (as in the simulations of McConnachie et al. 2009), explaining the increase SFR in M33 at that time. Furthermore, M32 could have been in the process of being stripped and both events could have contributed to the elevated star formation in M31.

Another possible explanation for the episode is a merger with a galaxy $\gtrsim$20% of the pre-merger mass of M31, such as in the merger simulations of Cox et al. (2008). The density profile of the stars formed, the recent decrease in the M31 star formation rate, and the stellar mass produced are all consistent with a relatively large merger scenario, as described below.

Like the Cox et al. (2008) simulations, the density profile of the enhancement is similar in shape to the overall surface density profile of M31, as shown in Figure 6. Furthermore, the simulations show a smooth decrease in global star formation rate down to levels $\lesssim$1 ${{M}_{\odot }}$ yr⁻¹ in the 2 Gyr following the merger as a result of gas depletion. Such a decrease is consistent with that recently measured by Lewis et al. (2015). The stellar density in the outer regions of merger simulations (e.g., Cox et al. 2008; Moreno et al. 2015) ∼3 Gyr after the merger show only relatively faint structures, qualitatively similar to those observed in the M31 halo (Ibata et al. 2001; Ferguson et al. 2002; McConnachie et al. 2009).

The amount of stellar mass produced by the event hints at a relatively large merger. If the matter that formed the stars in the burst came from the galaxy that merged, the baryonic Tully–Fisher relation (Zaritsky et al. 2014, e.g.,) yields a circular velocity of ∼140 km s⁻¹ for 2 × 10¹⁰ ${{M}_{\odot }}$ of baryonic mass, which is similar to the circular velocity of M33 (Corbelli & Salucci 2000). As some of the mass was in M31 already, this estimate is an approximate upper-limit on any infalling galaxy.

It is difficult to compare our SFHs of relatively low time resolution to the SFHs of the Cox et al. (2008) simulations directly. However, a 1:1 merger ratio in their simulations appears to produce more than the equivalent of a factor of 2 increase in star formation rate for 1.5 Gyr. Such an equal-mass merger appears to produce the equivalent of about a factor of $\gtrsim$10 increase for 1 Gyr. The 5.8:1 merger does not appear to enhance the star formation enough to produce the effect we see in the M31 data. The mean rate increase of a factor of 4–5 over about a Gyr we see in our data is limited by our time resolution. The actual burst could have had a much shorter duration, but produced enough stellar mass for us to detect it throughout the galaxy. Most consistent with merger simulations would be another galaxy $\gtrsim$20% of the pre-merger mass of M31. Such a total mass (∼2 × 10¹¹ ${{M}_{\odot }}$) is reasonable considering the baryonic mass estimate above of ∼2 × 10¹⁰ ${{M}_{\odot }}$.

The merger scenario also has problems. First, it is unclear if the M31 disk would have survived, even in its relatively quiescent state, if it had undergone such a merger. However, simulations do suggest that disks can survive large mergers if they contain high gas fractions (Hopkins et al. 2009). Second, it is puzzling why M33 would show a simultaneous global star formation episode. However, it is possible that the galaxy that merged with M31 passed by M33 (which may have been much closer to M31 at the time McConnachie et al. 2009) along the way, inducing star formation during the same epoch but not strongly disturbing its morphology. Finally, very recent simulations by Moreno et al. (2015) suggest that most of the star formation in mergers occurs <1 kpc from the galaxy center, which would not be consistent with this disk-wide burst in M31.

Whatever the explanation, these newly measured age distributions from the PHAT survey provide a new piece of evidence, along with the structured halo (Hammer et al. 2010, 2013; Sadoun et al. 2014), structured satellite distribution (Ibata et al. 2013, 2014a), peak in massive cluster ages (Fusi Pecci et al. 2005), and velocity dispersion of the disk (Dorman et al. 2015) that point to M31 having undergone a significant interaction event in the past few Gyr.

In conclusion, we have found new evidence that a global burst of star formation occurred in M31 2–4 Gyr ago over the entire M31 disk. The cause of this burst is far from clear. It could have been due to interactions with M32 and/or M33. It could have been a merger, or it could be some combination of these. The occurrence of the burst may help to explain the current, low-activity state of the M31 disk and its current global transformation toward the red sequence. Future simulations may be able to provide more detailed constraints on the physical cause of the burst.