The Hubble Space Telescope Survey of M31 Satellite Galaxies. II. The Star Formation Histories of Ultrafaint Dwarf Galaxies

Ultrafaint dwarf (UFD) galaxies around the Milky Way (MW) represent our strongest observational link between cosmic reionization and low-mass galaxy formation. Long-standing theoretical models posit that the ultraviolet (UV) background in the early Universe should prevent or suppress star formation in the lowest-mass dark matter halos (e.g., Efstathiou 1992; Bullock et al. 2000; Benson et al. 2002, 2003; Somerville 2002; Bovill & Ricotti 2009; Muñoz et al. 2009; Salvadori & Ferrara 2009; Busha et al. 2010; Tumlinson 2010; Simpson et al. 2013; Wheeler et al. 2015; Jeon et al. 2017). Stellar populations of these "fossil" galaxies were predicted to be largely ancient (z > 6), but color–magnitude diagram (CMD) analysis of classical Local Group dwarfs found virtually all of them to have substantial star formation at younger ages (e.g., Tolstoy et al. 1998; Gallart et al. 1999; Grebel & Gallagher 2004; Monelli et al. 2010b, 2010a; de Boer et al. 2012, 2014; Weisz et al. 2014a; Skillman et al. 2017; Savino et al. 2019). The discovery of UFDs in the early 2000s (e.g., Willman et al. 2005a, 2005b; Zucker et al. 2006a, 2006b; Belokurov et al. 2006, 2007; Sakamoto & Hasegawa 2006; Irwin et al. 2007; Walsh et al. 2007) and the precise determination of their star formation histories (SFHs) in the 2010s now provide the strongest evidence that reionization can stunt the formation of low-mass galaxies. Virtually all MW UFDs with SFHs measured from the oldest main-sequence turnoff (oMSTO) are consistent with a rapid quenching of their star formation at 6 ≲ z ≲ 8 (Brown et al. 2012, 2014; Okamoto et al. 2012; Weisz et al. 2014a, 2023; Gallart et al. 2021; Sacchi et al. 2021; Simon et al. 2021), in agreement with many predictions and observational constraints on the timing and duration of reionization (e.g., Fan et al. 2006; Ouchi et al. 2010; Pentericci et al. 2011; Ono et al. 2012; Robertson et al. 2015; Planck Collaboration et al. 2016).

Despite the uniformity of these SFHs, there remain lingering concerns about selection bias, as all of these UFDs share a common accretion history associated with the dark matter halo of the MW. While theoretical and observational efforts have been dedicated to quantifying environmental effects in the evolution of the MW's UFDs (e.g., Wetzel et al. 2015; Rodriguez Wimberly et al. 2019; Sacchi et al. 2021; Santistevan et al. 2023), the degree of synchronicity in formation and quenching among the dozens of known UFDs and the potential influence of local environment are poorly understood.

The cleanest way to alleviate these concerns is to measure the SFHs of UFDs outside the MW's halo. Ideally, the best test would use UFDs that are unambiguously isolated, mitigating any possible effects of environment. However, the first isolated UFD was discovered serendipitously and only within the past year (Sand et al. 2022). Instead, the only known UFDs that reside outside the MW halo, and that we can image down to the oMSTO, are satellites of M31 (e.g., Martin et al. 2016; Collins et al. 2022a; Martínez-Delgado et al. 2022; McQuinn et al. 2023). Measuring precise SFHs of these faint, distant systems (D ∼ 750 kpc) requires the excellent angular resolution and sensitivity of the Hubble Space Telescope (HST) and was a main science driver behind the 244 orbit HST Survey of M31 Satellites program awarded in Cycle 27 (HST-GO-15902, PI: Weisz).

In this paper, we present the SFHs of six UFDs that are in the satellite system of M31. These are the first SFHs of bona fide UFDs (i.e., M_V > − 7.5; Simon 2019) outside the MW's halo that have been measured from their oMSTOs. We summarize our data in Section 2, describe our methodology in Section 3, and discuss results in Section 4. Throughout the paper we assume a Λ cold dark matter (CDM) cosmological model based on Planck Collaboration et al. (2020).

We use photometric catalogs obtained as part of the HST Treasury Survey of the M31 Satellite System (Savino et al. 2022; D. Weisz et al. 2023, in preparation). This program obtained new F606W and F814W Advanced Camera for Surveys (ACS) and UVIS imaging of 23 satellites of M31 without previous deep HST imaging and was combined with archival HST data to provide oMSTO-depth imaging for 34 known M31 satellites.

The photometry for the UFDs considered in this paper is based on a combination of new F606W and F814W HST/ACS imaging from the HST Treasury program and shallower archival imaging in the same filters taken as part of program GO-13699 (PI:Martin; Martin et al. 2017). We select every galaxy targeted by the survey that has an absolute magnitude M_V > − 7.0, which leaves no ambiguity about their status as UFDs. We use absolute magnitudes from Savino et al. (2022), which are based on the PandAS structural parameters of Martin et al. (2016) and are recalibrated on the basis of updated RR Lyrae distances. This selection yields six targets: And xi (M_V = − 6.4 ± 0.4, 11 orbits of imaging), And xii (M_V = − 6.6 ± 0.5, 11 orbits), And xiii (M_V = − 6.8 ± 0.4, 9 orbits), And xx (M_V = − 6.4 ± 0.4, 11 orbits), And xxii (M_V = − 6.4 ± 0.4, 18 orbits), and And xxvi (${M}_{V}=-{6.0}_{-0.5}^{+0.7}$, 16 orbits).

Details of the photometric reduction and catalog construction are broadly described in Savino et al. (2022) and will be presented at length in the upcoming program survey paper (D. Weisz et al. 2023, in preparation). Here, we provide a brief summary. We reduced the ACS imaging using DOLPHOT (Dolphin 2000, 2016), a well-tested point-spread function photometry package that is commonly used in the analysis of HST imaging of nearby galaxies (e.g., Holtzman et al. 2006; Dalcanton et al. 2009, 2012; Monelli et al. 2010b; McQuinn et al. 2010; Radburn-Smith et al. 2011; Weisz et al. 2014b; Williams et al. 2014, 2021). We reduce the data using the same DOLPHOT setup recommended in PHAT (Williams et al. 2014), with the exception of the PSFPhotIT parameter, which we set to 2. This choice adds a second iteration to the photometric reduction, which refines the noise estimates and results in higher completeness at low signal-to-noise ratio (S/N).

From the resulting DOLPHOT catalogs, we select stars that are likely galaxy members using several criteria. First, we cull our photometric catalogs using DOLPHOT quality metrics. Namely, we select sources with:

• 1.  
  S/N ≥ 4
• 2.  
  Sharp² ≤ 0.2
• 3.  
  Crowd ≤ 0.75
• 4.  
  Round ≤ 3

in both F606W and F814W filters. This step eliminates spurious sources, including artifacts and extended background galaxies. These selection criteria are slightly different from the PHAT recommendations (Williams et al. 2014). This is because the lower stellar density in our satellite galaxies, compared to the disk of M31, allows us to apply more stringent thresholds to the values of Crowd (the amount of contamination from neighboring sources) and Round (the ellipticity of the source's light profile). We can therefore exclude contaminants more efficiently while retaining most stellar sources.

Second, we remove residual artifacts originating from the diffraction spikes of bright foreground stars. We do this by masking the region around all known Gaia stars with G < 18 that are located in or near our ACS fields. The masks are designed on a star-by-star basis to encompass the central saturated core and the visible extent of the diffraction spikes. Finally, we use the structural parameters of Martin et al. (2016) to select only sources that fall within the ellipse corresponding to 2 half-light radii (r_h ) from the photometric center of the galaxy. Our HST fields are close to the photometric center in each of our targets. In most cases, however, the ACS field of view is not sufficiently large to capture the full 2r_h ellipse. This means that the spatial selection function varies slightly among our targets. We do, however, fully account for this variable spatial footprint whenever relevant, e.g., to calculate stellar masses (see Section 3).

Figure 1 shows the CMDs for our six UFDs. The CMDs are deep, with F606W detection limits (S/N = 4) between 29.4 and 29.8. The S/N at oMSTO (M_F606W ∼ + 4.0) is ∼10. The stellar populations of these systems can be visually appreciated by the diversity of horizontal branch (HB) morphologies. We find that some of our galaxies have predominantly red HBs, confirming what was already reported from shallower HST photometry (Da Costa et al. 1996, 2000, 2002; Martin et al. 2017). This is particularly the case for And xiii and And xxii, which almost entirely lack a blue HB. And xx and And xxvi also have a high fraction of red HB stars. The red HB is typically associated with younger ages and/or higher metallicities (e.g., Gratton et al. 2010; Salaris et al. 2013; Savino et al. 2019). We do not include the HB in our SFH fits, and instead use the HBs as a secondary check on the results, which we discuss in Section 4.3.

Zoom In Zoom Out Reset image size

We model each UFD CMD using MATCH (Dolphin 2002), a commonly used software routine that recovers the system SFH by forward modeling the CMD with a combination of simple stellar population models. Details on the CMD-fitting methodology are provided in the original papers and in many other nearby galaxy studies (e.g., Weisz et al. 2014a; McQuinn et al. 2015). Here, we provide specific details used for modeling our sample.

We use a Kroupa initial mass function (Kroupa 2001), normalized between 0.08 M_⊙ and 120 M_⊙, and an unresolved binary fraction of 0.35, with secondary masses drawn from a uniform mass ratio distribution. The binary fraction in UFD galaxies is poorly known. However, because our CMDs do not extend significantly below the oMSTO, this parameter has very little impact. To verify this, we have repeated our analysis assuming binary fractions of 0.5 and 0.7, and obtained virtually indistinguishable results (more details are available in Appendix B). We assume homogeneous RR Lyrae based distances, anchored to Gaia eDR3, from Savino et al. (2022) and foreground extinction from Green et al. (2019). Our adopted distance and extinction values are listed in Table 1.

Table 1. Adopted Distance and Reddening Values for Our CMD Modeling, Literature Absolute Magnitudes, and Stellar Masses and Star Formation Timescales Inferred from Our Analysis of the Six UFDs

ID	(m − M)0	E(B − V)	MV	M*(z = 0)	M*(z = 5)	τ50	τ80	τ90	τ20–80	τ20–90
				(105 M⊙)	(105 M⊙)	(Gyr)	(Gyr)	(Gyr)	(Gyr)	(Gyr)
And xi	24.38 ± 0.07	0.09 ± 0.03	−6.4 ± 0.4	${1.4}_{-0.3}^{+0.3}$	${1.5}_{-0.8}^{+1.0}$	${13.2}_{-0.6}^{+0.9}$	${12.7}_{-0.1}^{+1.4}$	${9.8}_{-1.6}^{+2.8}$	${1.0}_{-1.0}^{+0.4}$	${1.2}_{-1.2}^{+4.4}$
And xii	24.28 ± 0.08	0.17 ± 0.03	−6.6 ± 0.5	${2.6}_{-0.9}^{+2.0}$	${2.6}_{-1.2}^{+3.0}$	${13.2}_{-0.6}^{+1.0}$	${12.6}_{-1.7}^{+1.5}$	${6.8}_{-0.6}^{+1.4}$	${1.1}_{-1.1}^{+1.8}$	${7.2}_{-2.1}^{+0.5}$
And xiii	24.57 ± 0.07	0.14 ± 0.03	−6.8 ± 0.4	${1.1}_{-0.2}^{+0.4}$	${0.1}_{-0.1}^{+0.2}$	${10.6}_{-0.6}^{+0.6}$	${10.1}_{-0.1}^{+1.1}$	${5.9}_{-0.4}^{+0.5}$	${0.9}_{-0.9}^{+0.2}$	${5.6}_{-1.4}^{+0.2}$
And xx	24.35 ± 0.08	0.08 ± 0.03	−6.4 ± 0.4	${1.0}_{-0.2}^{+0.2}$	${0.9}_{-0.5}^{+0.6}$	${13.1}_{-0.5}^{+1.0}$	${12.3}_{-1.7}^{+0.4}$	${9.1}_{-0.8}^{+1.6}$	${1.1}_{-1.1}^{+2.1}$	${5.1}_{-2.3}^{+0.5}$
And xxii	24.39 ± 0.07	0.10 ± 0.03	−6.4 ± 0.4	${1.1}_{-0.04}^{+0.5}$	${0.9}_{-0.2}^{+0.4}$	${13.0}_{-0.4}^{+1.1}$	${10.6}_{-0.6}^{+1.2}$	${8.1}_{-2.7}^{+2.1}$	${3.1}_{-1.6}^{+0.8}$	${3.6}_{-1.1}^{+4.8}$
And xxvi	24.48 ± 0.07	0.09 ± 0.03	$-{6.0}_{-0.5}^{+0.7}$	${1.8}_{-0.1}^{+1.3}$	${1.8}_{-0.7}^{+1.7}$	${13.1}_{-0.5}^{+1.0}$	${12.0}_{-1.0}^{+0.7}$	${5.7}_{-0.3}^{+1.0}$	${1.1}_{-1.1}^{+1.7}$	${7.3}_{-1.2}^{+1.2}$

Download table as:  ASCIITypeset image

We adopt the BaSTI scaled-solar stellar models (Hidalgo et al. 2018). While many galaxies in this mass range exhibit various degrees of α-enhancement (e.g., Simon 2019), this generally has a minor impact on the broadband filters, and scaled-solar models are used in many SFH studies (e.g., Monelli et al. 2010b; Skillman et al. 2017). The main difference is a zero-point offset in the recovered [Fe/H] values (Salaris et al. 1993; Cassisi et al. 2004). We show the effects of α-enhanced versus scaled-solar mixtures in Appendix B and find negligible differences in the recovered SFHs. The adequacy of BaSTI scaled-solar models is further illustrated in Appendix D, in which we show that the CMD of the metal-poor MW globular cluster M92 is well matched by the scaled-solar isochrone of the appropriate age and metallicity.

We use a grid of simple stellar populations over an age range of $7.50\lt {\mathrm{log}}_{10}(t/\mathrm{yr})\lt 10.15$ with a 0.05 dex resolution and metallicities that span −3.0 < [Fe/H] < 0, with a 0.1 dex resolution. We adopt a physically motivated prior on the age–metallicity relationship consistent with several past studies (e.g., Weisz et al. 2014a; Skillman et al. 2017). Specifically, we require the metallicity to increase monotonically with time, with a modest dispersion allowed at each age (0.15 dex). This helps mitigate the age–metallicity degeneracy at the oMSTO, which is affected by the modest S/N of the data, the limited temperature sensitivity of F606W-F814W, and the paucity of stars on the CMDs of these faint galaxies. This prior only requires metallicity to increase with time but puts no constraints on the metallicity values themselves. Additional details about the metallicity of our targets are provided in Section 4.3, and further tests on the effect of this assumption are detailed in Appendix B.

We model observational effects (photometric errors and incompleteness) for each galaxy using ∼5 × 10⁵ artificial star tests (ASTs). The ASTs are distributed uniformly in color and magnitude over the CMD and spatially distributed following a 2D exponential density profile that is taken from Martin et al. (2016). The ASTs are injected into the images and recovered with the same reduction setup described in Section 2.

We model the CMD from ∼2 mag above the brightest observed red giant branch stars down to the observed magnitude that corresponds to a 50% completeness level (28.3–28.6 in F606W and 27.6–27.9 in F814W). We use a CMD bin size of 0.05 × 0.1 mag in color and magnitude, respectively. We exclude the HB region from the fit due to uncertainties in the models of HB stars and to the availability of more reliable age and metallicity sensitive MSTO stars in our photometry (e.g., Gallart et al. 2005). In our models, we include MW foreground stars using the models from de Jong et al. (2010).

MATCH maximizes a Poisson likelihood function to find the most likely SFH that describes the observed CMD. We then calculate random uncertainties in the solution using the methodology described in Dolphin (2013), which is based on Hamiltonian Monte Carlo sampling of the solution parameter space. We also estimate the size of systematic uncertainties due to the stellar models, using the methodology of Dolphin (2012). This technique introduces perturbations in the luminosity and temperature of the reference stellar models, to simulate uncertainties in the stellar evolution parameters. We also explore the effects of stellar evolution models on our fits in Appendix B, and find they are consistent with our computed systematic uncertainties.

As an example for our modeling procedure, Figure 2 shows the observed CMD Hess diagram, the best-fit CMD model and the fit residuals for one of our targets (And xx). Aside from a minor discrepancy at the level of the HB (which was excluded from the fit; red box), there is good agreement between the observed and model CMD; no notable structure is visible in the residuals. The fits for the other five galaxies are of similar quality.

Zoom In Zoom Out Reset image size

We use the present-day and early Universe stellar masses as measured from CMD modeling as part of our analysis. To get the total stellar mass, we first integrate the SFH, which yields the birth mass of the stellar population in our HST field. The present-day stellar mass in our HST field is then derived by subtracting the amount of mass lost due to stellar evolution. We calculate this factor using the stellar population framework of Conroy et al. (2009), assuming a 13 Gyr stellar population with [Fe/H] = − 2.0 (see Section 4). The resulting mass loss is 42.5% of the stellar population birth mass. For old stellar populations, this value has very little dependence on the SFH. For stellar populations with 8 < t < 14 Gyr, and −3.0 < [Fe/H] < − 1.0, this factor changes by less than 2%. Then, we must account for the limited spatial coverage of the CMDs (due to the ACS field of view and the catalog spatial cuts). Under the assumption that the stellar mass-to-light ratio does not vary significantly across the galaxy, we use the morphological parameters from Martin et al. (2016) to estimate what fraction of the stellar mass falls outside our footprint. This fraction varies from 15% (And xx) to 55% (And xii). We then obtain the total present-day stellar mass for our UFDs.

The stellar mass uncertainties are calculated using a Monte Carlo approach. We obtain 5000 random realization of the stellar mass by sampling from the SFH uncertainties and the reported uncertainties in the Martin et al. (2016) structural parameters. We then use the 16th and 84th percentile of the stellar mass distribution to estimate our confidence interval.

We also estimate the stellar mass of the galaxies at z = 5 using a similar procedure. In this case, we use the total star formation in our oldest age bin to calculate the birth stellar mass. Because we have no constraints on when star formation occurred within our oldest time bin, stellar evolution mass loss is more uncertain. We include this effect in our Monte Carlo samples by drawing a stellar population birth time from a uniform $10.10\leqslant {\mathrm{log}}_{10}(t)\leqslant 10.15$ probability distribution and calculating the corresponding mass loss at z = 5. The median mass loss we obtain with this method is 34%. The resulting stellar masses of the galaxies at z = 0 and z = 5 are listed in Table 1. Because star formation causes stellar mass to increase over time and stellar evolution mass loss causes a decrease in stellar mass, the inferred stellar mass at z = 5 can be higher or lower than what is measured at z = 0, depending on the interplay between these two factors.

4.1. Lifetime SFHs

Figure 3 shows the lifetime cumulative SFHs for our six UFDs. We focus on the cumulative SFHs because they provide a robust statistical treatment of the uncertainties, which includes the covariance among different star formation bins. However, for completeness, we also provide the instantaneous star formation rates in Appendix A. Five of the UFDs in our sample (And xi, And xii, And xx, And xxii, and And xxvi) have a prominent episode of star formation in the earliest time bin (t > 12.6 Gyr), which drastically decreased in intensity by z ∼ 5. Taken at face value, the best-fit SFHs indicate that all five galaxies formed ≳50% of their stellar mass in the oldest time bin. This initial strong episode of star formation is followed by varying degrees of star formation at later times ranging from ∼15% (And xi) to ∼35% (And xxii) of the total stellar mass. The significance of this SFH tail, and its fidelity to the true evolution of our targets, will be further discussed in Section 4.3.

Zoom In Zoom Out Reset image size

And xiii exhibits a very different SFH. From the best-fit SFH, we see that only 10% of its stellar mass formed in the oldest time bin, followed by a ∼1 Gyr quiescent period. The galaxy then experienced a strong burst of star formation beginning 11 Gyr ago, in which it formed 75% of its total stellar mass in the span of 1 Gyr.

As we discuss below, the SFH of And xiii is unlike any other known UFD. Its unusual nature, and the observational challenges linked to this target (being the most distant system in our sample, one of the most extincted, and with relatively shallow photometry) motivated a number of tests for robustness (e.g., sensitivity to adopted stellar model, distance/extinction variations), all of which confirm this "young" UFD scenario. These tests are illustrated in the Appendices B, C, D, and E of this paper and reveal that, under virtually any assumption, no more than 40% of the total stellar mass of And xiii could have formed in the oldest time bin.

Figure 4 illustrates the features in And xiii's CMD that are linked to a younger age. We have taken the Hess diagram of And xiii's CMD, removed the effect of distance and reddening, and normalized its integral to 1, effectively deriving a probability density distribution of its stars on the CMD. We did the same for the Hess diagram of And xii, which has 80% of its star formation in the oldest time bin. We show the difference between the two density maps (left panel). Because the two photometric catalogs have very similar depth and completeness, differences in the density distribution are intrinsic to the stellar population properties.

Zoom In Zoom Out Reset image size

For comparison, we have also constructed a similar map (right panel) showing the difference between two stellar population models. The first one is an "early" star-forming model, that only has star formation in the oldest time bin ($10.10\lt \mathrm{log}(t)\lt 10.14$, with an average age of 13.5 Gyr), while the second is a "late" star-forming model, in which the star formation happens in our third time bin ($10.0\,\lt \mathrm{log}(t)\lt 10.05$, with an average age of 10.6 Gyr). The metallicities of the early and late models have been set to match the best fit for And xii ([M/H] ∼ − 1.75) and And xiii ([M/H] ∼ − 2.0), respectively.

As shown in Figure 4, the MSTO of And xiii is significantly brighter and bluer than that of And xii. This guides our fit to substantially younger ages. Although not included in the CMD fit, Figure 4 also highlights the much redder color of And xiii's HB compared to And xii, which is compatible with the younger SFH solution. In fact, both features (brighter MSTO and redder HB) are visible in the late-versus-early model comparison, although the models predict bluer HBs compared to observations, due to long-standing limitations in stellar mass-loss models (e.g., Gratton et al. 2010).

In principle, a higher fraction of ancient stars could exist in And xiii outside of our HST field of view. Stellar population gradients are known to be widespread in nearby classical dwarfs (e.g., Tolstoy et al. 2004; Battaglia et al. 2008; Vargas et al. 2014; Savino et al. 2019), with young metal-rich stars being more centrally concentrated than older, more metal-poor, stars. It is currently unclear whether these gradients exist in UFDs but, if they did, they would affect the fraction of young stars in our central HST field. However, the small size of And xiii (${r}_{h}\sim 0\buildrel{\,\prime}\over{.} 8$; Martin et al. 2016) means that our HST field contains ∼80% of And xiii's total light. Under the extreme assumption that every star outside our field of view was ancient, our SFH solution would still allow for no more than 30% of the total stellar mass to have formed before 12.6 Gyr ago.

4.2. Comparison with MW UFDs

To place our findings into a broader context, Figure 5 shows our best-fit SFHs, alongside similar measurements for the 13 MW UFD satellites that have published SFHs (Brown et al. 2014; Sacchi et al. 2021). The SFHs are color-coded according to the galaxy luminosity.

Zoom In Zoom Out Reset image size

Figure 5 shows two clear trends. First, 16 of the 19 galaxies in the combined sample formed ∼80% of their stellar mass by z ≥ 5. Two of the exceptions are Ursa Major I from Brown et al. (2014) and And xxii from this study, which formed 60% and 65% of their stellar mass by z = 5. The other exception is And xiii. Second, in general, the lowest-luminosity galaxies (M_V ≲ − 6) quench earlier than more-massive UFDs. This trend is clearly evident in the comparison with the Sacchi et al. (2021) sample, which is composed of particularly faint galaxies, but is also appreciable in the more-massive sample of Brown et al. (2014). Importantly, because the lowest-luminosity systems are typically MW UFDs and the higher-luminosity systems are M31 satellites, it is possible that this is not exclusively a stellar mass trend, but could be related to differences in the early environment of the two samples. We revisit this point later in the discussion.

Figure 6 shows the same comparison of Figure 5, color-coding the SFHs by stellar mass at z = 5. Because of the mostly ancient SFHs, the present-day luminosity roughly traces the high-redshift stellar mass, preserving the trend observed in Figure 5. However, And xiii stands out in this comparison. In fact, while having the highest present-day luminosity of our combined sample, And xiii has a stellar mass at z = 5 that is comparable to some of the MW UFDs. Nonetheless, while the latter quenched rapidly, And xiii halted star formation at a much later time (t ≲ 10 Gyr). We discuss the SFH of this galaxy in greater detail in Section 4.5.

Zoom In Zoom Out Reset image size

The trends of Figure 5 can be quantified more clearly through a comparison of characteristic star formation timescales. This is a common practice in the field (e.g., Weisz et al. 2014a, 2019; Skillman et al. 2017; Sacchi et al. 2021), which makes use of the quantity τ_X ; this is defined as the lookback time at which the galaxy reached X% of the total star stellar mass formed. For clarity, an illustration of these metrics is shown in Figure 5. Measurements of select τ_X for our sample are provided in Table 1.

Figure 7 compares values of τ_X for the same MW and M31 galaxies shown in Figure 5. It is again clear that, with the notable exception of And xiii, all MW and M31 UFDs formed at least 50% of their stellar mass by the end of cosmic reionization. Following previous studies (e.g., Skillman et al. 2017; Weisz et al. 2019; Sacchi et al. 2021) we can use τ₉₀ as a tracer of the quenching epoch (Figures 5(a), (b)). This metric reveals a clear difference between MW and M31 UFDs. MW UFDs quenched rapidly after reionization, whereas M31 UFDs sustained star formation until as late as z ∼ 1. These extended SFHs can be appreciated in Figure 5(b), which shows τ₂₀₋₉₀, defined as the time elapsed between 20% and 90% of the total star formation. While the MW UFDs have τ₂₀₋₉₀ of the order of 2 Gyr, M31's satellites tend to have significantly larger values ranging from 4–8 Gyr.

Zoom In Zoom Out Reset image size

In stark contrast to all other UFDs is And xiii. It has a younger stellar population with τ₅₀ = 10.6 ± 0.6 Gyr. Like other M31 UFDs in our sample, And xiii shows signs of having some residual star formation ongoing for several gigayears, after the major event of star formation (which happened between 10 and 11 Gyr ago).

4.3. Do M31's UFD Satellites Have Extended SFHs?

4.4. The Quenching of Local Group UFDs in the Context of Reionization

4.5. The Case of And xiii: The First Observed Episode of Reignition in an Ultrafaint Dwarf

4.6. Future Prospects

Using deep HST imaging that reaches the oMSTO, we derived lifetime SFHs for a sample of six UFD galaxies associated with the satellite system of M31. These are the first oMSTO-based SFHs for UFD galaxies outside of the MW halo. The main takeaways of our analysis are as follows:

• (i)  
  Five of the UFDs in our sample formed at least 50% of their stellar mass by z = 5 (12.6 Gyr ago), after which star formation activity drastically reduced in intensity. This is similar to other known UFDs and is compatible with the predicted effect of cosmic reionization. However, we find that star formation is not entirely quenched in our galaxies, and 10%–40% of their stellar mass formed at later times.
• (ii)  
  One UFD in our sample (And xiii) has a remarkable SFH, having formed only ∼10% of its stellar mass by z = 5 and having experienced a delayed period of strong star formation at 2 < z < 3. This UFD is the first of its kind and suggests that star formation reignition, previously observed in other more-massive galaxies, can also occur at lower stellar masses.
• (iii)  
  Combining SFHs for MW and M31 UFDs, we conclude that the effect of reionization is not homogeneous among UFDs, with more-massive galaxies being able to sustain low levels of prolonged star formation compared to low-mass UFDs, which are quenched more rapidly. However, due to selection biases, the environmental dependence of this trend cannot be ruled out.
• (iv)  
  Future observational efforts, including large-area surveys, and dedicated follow-ups with HST and JWST, will provide an extended observational baseline to address some of the remaining open questions in UFD evolution. These include a better understanding of the role of cosmic environment, establishing the degree of SFH diversity in galaxies at the massive end of the UFD spectrum, and bracketing the stellar/halo mass range where the transition between ancient reionization fossils and classical dwarf galaxies occurs.

The authors thank the anonymous referee for providing a constructive report. A.S. thanks E. Sacchi for graciously sharing the star formation histories of seven Milky Way satellites. Support for this work was provided by NASA through grants GO-13768, GO-15476, GO-15902, AR-16159, and GO-16273 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. A.W. acknowledges support from: NSF via CAREER award AST-2045928 and grant AST-2107772; NASA ATP grant 80NSSC20K0513; HST grants AR-15809, GO-15902, GO-16273 from STScI. M.C.C. acknowledges support though NSF grant AST-1815475. M.B.K. acknowledges support from NSF CAREER award AST-1752913, NSF grants AST-1910346 and AST-2108962, NASA grant NNX17AG29G, and HST-AR-15006, HST-AR-15809, HST-GO-15658, HST-GO-15901, HST-GO-15902, HST-AR-16159, and HST-GO-16226 from STScI. E.N.K. acknowledges support from NSF CAREER grant AST-2233781. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

Facility: HST - Hubble Space Telescope satellite (ACS). All of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via doi:10.17909/ftep-8k64.

Software: This research made use of routines and modules from the following software packages: Astropy (Astropy Collaboration et al. 2013), DOLPHOT (Dolphin 2016), IPython (Perez & Granger 2007), MATCH (Dolphin 2002), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), Pandas (McKinney 2010), and SciPy (Virtanen et al. 2020).

Here we provide the instantaneous star formation rates, also referred to as differential SFHs, derived from our fiducial solutions. These are illustrated in Figure 8, alongside the corresponding cumulative SFHs. Differential SFHs offer an intuitive way to visualize the star formation activity of the galaxy at any point in time. However, it is important to exercise caution when interpreting the uncertainties associated with these measurements. The uncertainties in the histograms, in fact, come from the variance of the star formation rate within each time bin. On the other hand, due to the details of the CMD-fitting method (see, e.g., Dolphin 2002), significant covariance exists among adjacent star formation bins. This constrains the range of plausible solutions that accurately fit the CMD.

Zoom In Zoom Out Reset image size

For instance, in Figure 8, most of the late star formation bins have error bars compatible with zero star formation (see, e.g., And xx). However, the star formation rate cannot be zero simultaneously in all bins, as that would significantly impact the star counts within the CMD. To address this limitation, cumulative SFHs are conventionally preferred, because they are constructed to incorporate the covariance between star formation rates. In the case of And xx, the statistical uncertainties on the cumulative SFH (yellow shading), unequivocally reject a complete absence of star formation after the oldest time bin.

One of the primary sources of systematic uncertainties in the CMD modeling of external galaxies lies in the choice of stellar evolution models (e.g., Gallart et al. 2005). In the main text, we have estimated the magnitude of these effects by using the methodology of Dolphin (2012), which generates a large set of mock stellar models by applying systematic shifts in luminosity and effective temperature to the adopted stellar evolutionary tracks. Here, we provide a different estimation of the stellar evolution systematics, by rerunning our analysis with a range of different stellar evolution models. In addition to our fiducial stellar model grid (BaSTI; Hidalgo et al. 2018), we have modeled the CMDs of our targets using the PARSEC (Bressan et al. 2012; Chen et al. 2014), MIST (Choi et al. 2016; Dotter 2016), and Victoria–Regina (VandenBerg et al. 2014) stellar models. Furthermore, while our fiducial analysis has been done using a scaled-solar abundance mixture, here we assess the effect of the assumed α-enhancement by including the [α/Fe] = 0.4 models of the BaSTI (Pietrinferni et al. 2021) and Victoria–Regina grids.

The resulting SFHs are shown in Figure 9. First, we note that the adopted α-enhancement has little effect on the shape of the SFH. The scaled-solar and α-enhanced solutions are virtually indistinguishable in most cases. In this regard, the major difference that the α-enhancement makes is to introduce a zero-point shift in the [Fe/H] scale of our solutions. This is a well-known effect from stellar evolution theory (e.g., Salaris et al. 1993). Second, the range spanned by the different SFH solutions is well-captured by the systematic confidence interval from the Dolphin (2012). The latter uncertainty interval tends to be slightly larger compared with the test presented here. But, as discussed in Dolphin (2012), the systematics are designed to be conservative given that there are only a limited set of stellar libraries available for use; they are unlikely to include the full range of reasonable stellar evolution uncertainties. Finally, the qualitative conclusions drawn from our SFHs are broadly unaffected by the choice of stellar models. With the exception of And xiii, all of our galaxies have a predominantly ancient population, regardless of stellar model choice. However, all stellar libraries suggests that these galaxies experienced significant star formation (10%–40% of total stellar mass) after reionization, with the only exception of And xi, which appears purely ancient.

Zoom In Zoom Out Reset image size

For And xiii, three out of the four stellar libraries used here confirm that the galaxy had little to no star formation at ancient times. The BaSTI and PARSEC models date the start of significant star formation at around z ≃ 3, while the Victoria–Regina models put it around z ≃ 4. Only the MIST stellar library allows for the presence of ancient stars in this galaxy and, even then, at least 30% of the stellar mass must have formed at a later time. It must be noted that, while the MIST library estimates an older age for And xiii, compared to the other models, the same happens for the other five UFDs. For these galaxies, in the absence of our cosmological limit of 14 Gyr on the grid edge, the best-fit SFH from MIST would peak at ages much older than the age of the Universe. This means that, even with the MIST models, And xiii remains a younger galaxy compared to the rest of the UFD sample.

The values of τ₅₀, τ₈₀, and τ₉₀ measured from the different solutions of Figure 9 are listed in Table 2. The measured values of τ₅₀ are in very good agreement, with an average standard deviation of ∼200 Myr across the different stellar evolution libraries. The values of τ₈₀ and τ₉₀ have higher standard deviations of ∼700 Myr and ∼1.6 Gyr, respectively. These characteristic spreads among the models are comparable to the uncertainties listed in Table 1. Overall, the result of these tests is that, while the choice of stellar models can affect the details of the measured SFHs by a modest degree, the key SFH features identified in the paper are robust, regardless of stellar library adopted.

As noted in Section 3, we have also tested the effect of adopting different binary fractions, by rerunning our fits with a binary fraction of 0.5 and 0.7. These solutions are illustrated in Figure 10, along with our fiducial SFHs, which use a binary fraction of 0.35. The result of this test is that the choice of binary fraction has an almost negligible effect on the recovered star formation history, with only And xxii showing a modest difference in the amount of star formation in the oldest bin.

Zoom In Zoom Out Reset image size

An additional test on the stellar population models regards the prior on the age–metallicity relation (AMR). Our fiducial SFHs are based on the assumption that the AMR monotonically increases with time. This is a common practice (e.g., Weisz et al. 2011, 2014a; McQuinn et al. 2015) that results in physically reasonable solutions even at low S/N on the oMSTO. However, given the results of this paper, especially on And xiii, a monotonically increasing AMR might not necessarily be an appropriate assumption. For instance, the late star formation we measure in some of our targets could have been triggered by the accretion of low-metallicity gas. We have therefore calculated a set of alternative SFH solutions, in which the AMR is left free to assume any form.

The results of this test are shown in Figure 11. Removing the constraints on the AMR has, for the most part, a modest effect. Most notably, And xiii remains a purely "young" galaxy, while And xi maintains its predominantly ancient SFH. The effect on the other targets is more pronounced, especially in And xii. However, changing our assumptions on the AMR does not impact the main conclusion of the paper, i.e., the presence of post-reionization star formation in some UFDs of our sample. On the contrary, the alternative SFHs of Figure 11 have a higher amount of late star formation, compared to our fiducial solutions. However, we advise strong caution in interpreting these results. By removing the requirement of a monotonic AMR, we have dramatically increased the degrees of freedom in our solutions, which can result in an unphysical metallicity evolution. In the case of Figure 11, the new solutions have strongly fluctuating AMRs, with values of [M/H] approaching solar values and varying by more than 1 dex between adjacent time bins. Nevertheless, a safe conclusion from this experiment is that our assumptions on the AMR do not impact the presence of late star formation in our sample, especially with respect to the case of And xiii.

Zoom In Zoom Out Reset image size

Two assumptions made in SFH measurements are the distance of the observed galaxy and the amount of foreground dust extinction. These will affect the magnitude and color of the observed stars and, potentially the SFH inferred by CMD fitting. To assess the effect of these parameters on our SFH measurements, we have repeated the analysis varying the assumed distance and reddening values, in accordance with the reported uncertainties in the literature. We have adopted distance uncertainties from Savino et al. (2022). For the six dwarf galaxies of our sample, these are of the order of 0.07–0.08 mag. Reddening uncertainties were taken from the dust maps of Green et al. (2019) and amount to 0.03 mag in E(B − V). We have then calculated a set of alternative SFHs, in a 3 × 3 grid, by applying [−1σ, 0, +1σ] shifts in distance and [−1σ, 0, +1σ] shifts in reddening to our CMD models (the central point of this grid is the trivial case where no shifts are applied and corresponds to our fiducial SFH solution.).

The results of these tests are illustrated in Figure 12 as blue lines and compared with the fiducial SFH used in the paper (black). The effect of distance and reddening uncertainties is comparable in size with other sources of systematic uncertainties included in our error budget (gray region in Figure 12). They modestly change the details of the SFHs but overall preserve the results presented in the paper. Most of the solutions for And xi, And xii, And xx, And xxii, and And xxvi have a main burst of star formation at the oldest ages, followed by more recent star formation accounting for 10%–40% of the total stellar mass. Conversely, And xiii maintains a significantly delayed star formation compared to the other dwarfs.

Zoom In Zoom Out Reset image size

We note from Figure 12 that certain combinations of distance and reddening can result in substantial changes in the inferred SFH. For the cases of And xii, And xx, and And xxvi this manifests as a significantly younger solution compared to the other setups. These solutions correspond to the extreme points of our parameter space, where both distance and reddening are shifted by 1σ. Under the assumption of Gaussian uncertainties, the chance of this happening for a single galaxy is ∼10% and much lower for the entire set of galaxies.

A clear difference between our observations of M31 satellites and the data used to derive the SFH of MW UFDs is the photometric depth. Due to the distance of M31, our CMDs are shallower and of lower S/N at the MSTO than the exquisite CMDs that have been obtained within the halo of the MW. To test how this difference may affect SFH recovery, we have taken two well-studied nearby stellar systems and created mock observations as if they were at the distance of M31.

The results of this test are shown in Figure 13. We begin with deep HST observations of the globular cluster M92 (GO-9453, PI: Brown; GO-10775, PI: Sarajedini) and of the UFD Eridanus II (GO-14234, PI: Simon). These data have S/N values at the oMSTO of ∼600 (M92) and ∼40 (Eridanus II), which are both much larger than those of our M31 photometry (S/N ∼ 10). For this test, we first derived the CMD-based SFHs using the original photometry of M92 and Eridanus II. This is done using the same methodology outlined in the paper. We assume (m − M)₀ = 14.59 and E(B − V) = 0.02 for M92 (Harris 2010), and (m − M)₀ = 22.84 and E(B − V) = 0.01 for Eridanus II (Martínez-Vázquez et al. 2021). The resulting SFHs are the red lines in the right panel of Figure 13 and will be assumed as our ground truth. We also use the well-known properties of M92 to further validate our overall methodology. The SFH we recover for M92 essentially corresponds to a simple stellar population with age and metallicity in excellent agreement with literature values (e.g., Dotter et al. 2007; Harris 2010; VandenBerg et al. 2013; the difference between our best-fit [Fe/H] of −2.0 and the measured abundance of −2.35 is completely due to our use of scaled-solar models to fit an α-enhanced population). The corresponding best-fit isochrone is also shown, in Figure 13, on the M92 CMD to illustrate the fidelity of the BaSTI stellar models in the F606W/F814W photometric bands.

Zoom In Zoom Out Reset image size

We then shifted the photometry of these two targets to simulate the same distance and reddening of one of our galaxies (And xiii), i.e., (m − M)₀ = 24.45 and E(B − V) = 0.14. We then perturbed the measurements using the photometric uncertainties and the completeness estimated from the artificial star tests of And xiii. The resulting CMDs (left panel of Figure 13) simulate how the M92 and Eridanus II CMDs would appear if they were satellites of M31 observed by our Treasury program. We repeat the CMD modeling for these noisy CMDs and find that, other than the larger uncertainties, the resulting SFHs (black lines in the right panel of Figure 13) are virtually indistinguishable from the ones obtained from the deep CMDs. The main difference is the increase in uncertainty, most notably for Eridanus II, owing to the lower S/N at the MSTO. This test therefore rules out any appreciable systematic in the SFHs caused by the limited S/N of our data at the MSTO.

We checked the potential impact of the choice of CMD-fitting software on our derived SFHs by comparing our solutions to those obtained with an independent methodology. We use the method described in Cole et al. (2007, 2014) to re-derive the SFHs of And xi and And xiii, which have the oldest and youngest stellar populations in our sample, respectively. We use identical photometric samples, stellar libraries, distance, and reddening as used in the main analysis to isolate the effect of the fitting methodology on the results. We summarize the setup of the "Cole" method below.

The BaSTI isochrones at each metallicity from Z = 0.0001–0.001 are binned into uniformly logarithmically spaced age groups from 9.70 $\leqslant \mathrm{log}(t)\leqslant$ 10.13 before being weighted by the initial mass function. Stars in each age bin are allowed to have any metallicity in the library, but the oldest bins ($\mathrm{log}(t)\geqslant 10.0$) are prevented from having Z > 0.0006. Each age bin is initially taken to be 0.10 dex wide, but may grow or shrink during the fitting depending on the noise level in the fit. We take the initial mass function from Chabrier (2003), which is log-normal within 0.08 ≤ M/M_☉ ≤ 0.8. Only 35% of the stars are assumed to be true singles and the rest are binaries; distinguishing between wide and close binaries, three-quarters of the binaries have the secondary mass drawn from the Chabrier initial mass function with the remainder drawn from a flat (uniform) initial mass function. The weighted isochrones are shifted by the adopted distance and reddening, convolved with the color and magnitude shifts from the artificial star tests, and weighted by the incompleteness fraction. The CMD is divided into bins that are 0.05 by 0.10 mag in size, and the SFH and AMR are derived by maximizing the Poisson likelihood for star counts in each bin that contains at least one star (excluding the same region of the HB and red clump that was excluded in the MATCH fits).

Figure 14 shows the results of the CMD fit, compared to the fiducial solutions, for And xi and And xiii. The Cole solution for And xi confirms that the majority of stars in this galaxy are compatible with the oldest ages allowed by the model grid. The differences between the Cole and MATCH solutions are entirely attributable to the different extension of the grid to old ages. The star formation in And xiii is distributed on a slightly longer episode in the Cole solution, compared to the narrower burst obtained with MATCH. Overall, the difference between the two And xiii solutions is mostly contained within the systematic uncertainty envelope (shaded gray region) and confirms the broad consistency with MATCH in previous studies of isolated dwarf galaxies (e.g., Monelli et al. 2010b, 2010a; Skillman et al. 2014). In spite of these subtle differences, the comparison between And xi and And xiii is unchanged. The Cole SFH of And xi reveals a primarily ancient stellar population, while And xiii shows prominent star formation well beyond reionization, with only ∼15% of the total star formation happening at t < 12.3 Gyr (z < 4). This test supports the robustness of our results against the choice of CMD-fitting methodology.

Zoom In Zoom Out Reset image size