The Hubble Space Telescope Survey of M31 Satellite Galaxies. III.
Calibrating the Horizontal Branch as an Age Indicator for Nearby Galaxies

We use HST/ACS F475W, F606W, and F814W photometric catalogs uniformly produced by the HST M31 Treasury Survey of M31 Satellites (GO-15902; PI: Weisz). The photometry used in this paper comes from new observations, taken as part of GO-15902, as well as additional archival data of the M31 system (GO-13028/13739, PI: Skillman; GO-13699, PI: Martin). The photometric reduction process with DOLPHOT (Dolphin, 2000) and subsequent catalog construction is summarized in Savino et al. (2022, 2023) and will be fully detailed in Weisz et al. (in prep). Here, we focus on highly complete and high signal-to-noise ratio (SNR) stars around the HB, which are robustly measured regardless of the specific choices in the reduction set-up. We use this M31 sample, as opposed to the Milky Way (MW), because of its uniformity of data quality, areal coverage, and SFH and relatively small degree of MW foreground contamination relative to the angularly larger MW satellites. Similarly, the chemical abundance peculiarities in MW globular clusters appear to influence the morphology of the HB in ways that are different from galaxy field populations, limiting the utility of globular clusters for our purposes (Gratton et al., 2011; Dalessandro et al., 2011; Bastian & Lardo, 2018, and references therein).

From the survey’s full sample of 36 dwarf galaxies (Savino et al., 2022), we exclude M32, NGC147, NGC185, NGC205, and And XIX as their large angular extent, off-center position of the HST pointings, and our incomplete knowledge of population gradients in these galaxies (e.g., Rose et al., 2005; Ho et al., 2015; Taibi et al., 2022) make it challenging to estimate representative metallicities that are matched to our fields (cf. § 2.3). Moreover, due to crowding, M32 and NGC205 do not have sufficiently deep photometry to ensure reliable stellar population ages (i.e., their CMDs do not reach the oMSTO).

We also exclude IC1613, the Pisces dwarf (LGS3) and the Pegasus dwarf irregular, as these galaxies have had recent star formation, and contamination from bright main sequence (MS) stars makes a clean HB selection more challenging for our calibration purposes. In fact, we use the Pisces and Pegasus dwarf galaxies in § 4.2, to validate how our methodology performs in the presence of main-sequence contamination of the HB. We also exclude And IX, which suffers from significant contamination from M31 itself, making it challenging to isolate genuine member HB stars.

Figures 1 and 2 show the CMDs of our final sample, focused on the HB region. Overall we have 21 galaxies with F606W/F814W data and 6 galaxies with F475W/F814W data. The photometry has been transformed to the absolute magnitude space using distances from Savino et al. (2022) and foreground extinction from Green et al. (2019). All the analysis in this paper will be carried out on absolute magnitudes. The typical SNR at the HB is between 60 and 150, and the completeness, as estimated from artificial star tests, is $\sim 100$%.

The stellar population ages used in this work are based on a set of homogeneous SFHs measured as part of the M31 Treasury survey (Savino et al., in prep; HST-GO-15902, PI D. Weisz). The general process of measuring these SFHs, along with the SFHs of the ultra-faint M31 satellites, are detailed in Savino et al. (2023).

In short, we derived the SFHs with the widely-used MATCH software package (Dolphin, 2002, 2012, 2013). MATCH models the density distribution of stars in the CMD, down to the oMSTO in the case of all our data. We used the same distance and extinction values used in this paper (Green et al., 2019; Savino et al., 2022). The CMDs were modeled using the BaSTI stellar model library (Hidalgo et al., 2018), covering an age range $10.15<\log{t}<7$, with bins of size $0.05$ dex, and a metallicity range of $-3.0<[Fe/H]<0.5$, with bins of size $0.1$ dex. The HB was explicitly excluded from the CMD fit and the age information we derive is primarily from the MSTO. The uncertainties on the SFH are then calculated using the methodology described in Dolphin (2012, 2013). Random uncertainties are based on Hamiltonian Monte Carlo sampling of the solution parameter space, while systematic uncertainties are estimated by introducing systematic perturbations in the shapes of the stellar models.

For each galaxy, we use the MATCH-based SFH to derive the mass-weighted age of the stellar population, $t^{*}$, defined as:

where $t_{n}$ is the median age in the $n^{th}$ star formation bin and $M^{*}_{n}$ is the total stellar mass formed in that bin. As the HB is only composed of low-mass stars, we only average over the star formation bins older than 6 Gyr, which corresponds to stellar masses on the HB of roughly $1M_{\odot}$. As the precise age upper limit probed by our analysis is somewhat uncertain, we experimented with other age cuts and found that 6 Gyr provides the strongest correlation with HB morphology. The values of $t^{*}$ obtained in this way are listed in Tab. 1. We calculate random and systematic uncertainties on $t^{*}$ by sampling from the respective SFH uncertainties.

As our analysis is based on ages from a homogeneous dataset and methodology, we only use the random uncertainties in our models. This is sufficient to achieve a robust relative age calibration of the HB morphology and explore galaxy-to-galaxy differences. However, the effect of systematic uncertainties needs to be taken into account when placing our HB-based ages in a broader context. Systematic uncertainties on the value of $t^{*}$, inferred from the MSTO, are listed in Tab. 1. While the precise value of this systematic term changes from galaxy to galaxy, the HB-based age scale we present in this paper is derived from the full sample of galaxies. For this reason, we adopt the median value of the systematic uncertainty distribution, which is 0.86 Gyr, as our systematic uncertainty on the HB-based ages.

At present, only a modest subset of M31 satellites, typically the brightest systems, have published spectroscopic metallicities of a large number of resolved RGB stars (Vargas et al., 2014; Ho et al., 2015; Kirby et al., 2020). We thus cannot rely on direct metallicity information for all galaxies in our sample. Instead, in order to maintain uniformity across the sample, we adopt a mean metallicity value for the stellar populations in our fields using the well-established luminosity-metallicity (LZ) scaling relation from Kirby et al. (2013b):

which has been obtained from spectroscopic studies of Milky Way dwarf spheroidal galaxies and that has been shown to be adequate for M31 satellites as well (Ho et al., 2015; Kirby et al., 2020). To calculate the metallicities we use stellar luminosities from Savino et al. (2022). The uncertainty on $\rm\langle[Fe/H]\rangle$ is calculated by propagating the uncertainties on the coefficients of eq. 2 and those on the luminosities, and adding the resulting term in quadrature with the LZ scatter of 0.17 dex reported in Kirby et al. (2013b). The values of $\rm\langle[Fe/H]\rangle$ obtained in this way are listed in Tab. 1. In Appendix A, we explore the effect of adopting LZ-based metallicities, instead of spectroscopic values, and find differences compatible with the uncertainties in our model.

The value of $\rm\langle[Fe/H]\rangle$ obtained through Eq. 2 is a mean metallicity of the galaxy stellar population (as measured from RGB stars, in the original study of Kirby et al., 2013b); it is not guaranteed to be representative of the HB population (e.g., Nagarajan et al., 2022; Savino et al., 2022), because HB and RGB stars do not come from identical stellar sub-populations and because of evolutionary timescale differences. However, direct metallicities of HB stars are impossible to obtain in any extragalactic system except the closest ones (e.g. Clementini et al., 2005). The next best approach is to rely on RGB-based tracers in both our calibration sample and in future applications to distant galaxies. This potential source of uncertainty is treated as a nuisance parameter in our models, as further discussed in § 3.2.

The morphology of the HB can be parametrized in a number of ways (e.g., Mackey & van den Bergh, 2005; Dotter et al., 2010; Milone et al., 2014). Here we adopt a commonly used (e.g, Da Costa et al., 1996, 2000, 2002; Martin et al., 2017) population ratio $\eta$, defined as:

where $n_{BHB}$ and $n_{RHB}$ are the number of blue and red HB stars, respectively. The selection of BHB and RHB stars is based on the Martin et al. (2017) definition, which makes use of CMD boxes. However, we have refined our box locations to accommodate our improved photometric precision, the use of multiple photometric bands, and to explicltly address contamination from the RGB and, if present, the MS. We now provide a detailed description of our HB member selection.

The left panel of Fig. 5 shows the mass-weighted age, $t^{*}$, of the old ($t>6$ Gyr) stellar population in our galaxy sample, as a function of $\eta$ and $\rm\langle[Fe/H]\rangle$ value. Consistent with expectations from stellar models and previous observational work (e.g., Sarajedini et al., 1995; Dotter et al., 2010; Salaris et al., 2013), we observe a general trend that galaxies with a red HB morphology (low $\eta$ values) tend to have younger stellar populations, while galaxies with blue HB morphologies (high $\eta$ values) are on average older. There is also indication that, at fixed $\eta$, galaxies of higher metallicity might be comparatively younger although the trend is less defined.

A clear takeaway from Fig. 5 is that galaxies with particularly blue HB morphologies ($\log(\eta)>-0.5$), all have predominantly old and metal-poor stellar populations, having $t^{*}>12$ Gyr and $\rm\langle[Fe/H]\rangle<-1.9$. This is in accordance with expectations from stellar evolution theory (e.g., Cassisi et al., 2013), as hot temperatures on the HB can only be reached by low-mass metal-poor models. As the HB morphology becomes redder, the scatter of the relation increases. This is the manifestation of the dependence of the HB temperature on age and metallicity (both mean and dispersion). Similar HB morphologies, in fact, can be obtained from different combinations of age and chemical composition.

Given this degeneracy, we model the relationship between $\eta$ and age (Fig. 5) using a bivariate function:

where $log(t^{*})$, $log(\eta)$ and $\rm\langle[Fe/H]\rangle$ are assumed to follow a linear relation with an intrinsic scatter of variance $V$. We fit this model using the affine invariant ensemble Markov chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al., 2013). The free parameters of this model are $a$, $b$, $c$, as well as $ln(V)$. We adopt uniform priors listed in Tab. 2.

We sample the posterior distribution using 32 walkers, setting the burn-in to first 100 steps, and defining the convergence length of the MCMC chain as 50 times the autocorrelation length. As both $log(\eta)$ and $log(t^{*})$ have asymmetric uncertainties, we calculate our model likelihood under the assumption that these quantities follow a split-normal distribution. Figure 6 shows the corner plot of our MCMC runs. Our posterior distributions are well-constrained and unimodal. Adopting the 50th percentiles as our fiducial parameters, and the 16th and 84th percentiles as our confidence interval, the fiducial relation is:

with intrinsic scatter $ln(V)=-7.671_{-0.572}^{+0.447}$, which corresponds to an rms of $0.022$ in $log(t^{*})$. For a population of 11 Gyr, this corresponds to a $0.55$ Gyr scatter. Our fiducial model is shown as a black line in the right panel of Fig. 5. Using eq. 7, the value of $log(t^{*})$ can be determined from the value of $\eta$ and $\rm\langle[Fe/H]\rangle$, with uncertainty:

where $\mathbf{C}$ is the covariance matrix

which is derived from our MCMC chains, and $\mathbf{J}$ is the Jacobian matrix of Eq. 7.

We comment on a few aspects of Eq. 7. The first regards the equation coefficients. The coefficient of the $\eta$ term is positive, meaning that bluer HB morphologies correspond to older stellar population ages. This is evident from Fig. 5 and is expected from stellar evolution. The coefficient of the metallicity term, on the other hand, is negative, meaning that, at higher metallicities, a given HB morphology traces a younger stellar population age. This seems at odds with stellar evolution expectations, as higher metallicities should require smaller masses (i.e., older ages) to achieve a given temperature on the HB. The interpretation of the $\rm\langle[Fe/H]\rangle$ term is more nuanced for a few reasons.

First, the metallicity obtained through Eq. 2 represents the stellar metallicity averaged across all stellar populations capable of producing an RGB. This is not guaranteed to be representative of the metallicity on the HB (which, given our selection, is produced only by stars older than $\sim 6$ Gyr), especially in galaxies with extended SFHs. In fact, the value of $\rm\langle[Fe/H]\rangle$, as obtained from Eq. 2, is essentially a proxy for stellar luminosity, rather than a direct probe of the metallicity of the old stars.

Second, and more importantly, ages and metallicities of stars in a galaxy are not truly independent variables. On the contrary, well defined age-metallicity relations (AMRs) are known to exist in nearby galaxies, as a result of chemical self-enrichment (e.g., de Boer et al., 2012a, b). These correlations exist within individual galaxies but also across the galaxy sample, with higher metallicity (i.e., higher luminosity) systems having on average more extended star formation histories and, therefore, younger mean ages. This is observed, for instance, in satellites of the MW (e.g., Weisz et al., 2015) and holds true for our M31 sample (Savino et al., in prep.). The negative coefficient of eq. 7 is likely capturing this correlation between galaxy mean age and metallicity. At fixed metallicity, the morphology of the HB then allows to quantify differences in mean stellar population age.

As this correlation between mean age and metallicity is encoded in our model, it follows that our calibration works best on galaxies that have broadly compatible AMRs to those of the M31 satellites. This could be a possible source of additional scatter when comparing galaxies across very different hosts. However, from eq. 7, this effect only becomes larger than the intrinsic scatter term if the galaxy AMR is offset by more than approximately 0.3 dex from those of our calibration sample.

Another interesting result in Eq. 7 regards the size of the intrinsic scatter term. The presence of a non-zero intrinsic scatter in our relation is not surprising, and this is likely arising from several sources of uncertainties in our methodology. These include the previously discussed discrepancy that could exist between $\rm\langle[Fe/H]\rangle$ and the HB metallicity. Another source of scatter may be that all the variables of our model ($t^{*}$, $\eta$, and $\rm\langle[Fe/H]\rangle$) describe global properties of the stellar population (being integrals of the age, HB color, and metallicity distributions, respectively). While stellar evolution suggests that a strong correlation should exist between age, metallicity, and HB color on a star-by-star basis, the correlation between global stellar population measurements might be somewhat reduced. For instance, two galaxies may have distinct stellar age distributions that result in similar values of $t^{*}$, while their HB morphologies might not necessarily lead to comparable values of $\eta$. Finally, we are assuming that the other variables affecting the HB (e.g., helium abundance or RGB mass loss) do not vary significantly across the galaxy sample, whereas some variation might be present.

Regardless, all these effects only introduce modest scatter in our fit. As shown in Fig. 5, the inclusion of a 0.022 dex $rms$ term in our model (0.5 to 0.7 Gyr, depending on the $t^{*}$) is sufficient to capture these additional uncertainties for virtually all galaxies (with a few exceptions discussed below). As demonstrated in § 4.4 and § 4.5, this intrinsic scatter is our dominant source of uncertainty on a large range of galaxy luminosities and photometric depths. The resulting age uncertainty is comparable in size to the typical systematics that affect old stellar population ages (e.g., Tab. 3), meaning that our HB calibration is robust enough to be used as an effective archaeological tracer.

The right panel of Fig. 5 demonstrates that our model is able to capture the relation between age, metallicity, and HB morphology, for most galaxies in our sample, within measurement uncertainties and the intrinsic scatter term. However, galaxies And VI, And XIII, and And XVI are clearly outliers relative to our model and most of the other systems. These three galaxies appear to be significantly younger than what we predict based on their HB morphology and metallicity.

Interestingly, these three galaxies also have unusual SFHs compared to the rest of the sample. As demonstrated in Savino et al. (2023), And XIII has experienced unusually low levels of early star formation, compared to other M31 satellites of similar luminosity, and ignited prominent star formation much later, at $z\lesssim 3$. A similarly delayed SFH has been shown for And XVI (Weisz et al., 2014b; Monelli et al., 2016; Skillman et al., 2017) and, from a recently measured SFH (Savino et al., in prep.), And VI also shows a strong burst of star formation at $z\lesssim 3$.

In the context where our model described by eq. 7 is capturing the correlation between $t^{*}$ and $\langle[Fe/H]\rangle$ introduced by the AMR, it is plausible that the late dominant star formation episodes for these three galaxies resulted in a similarly unusual enrichment history. For instance, the ignition of late prominent star formation could have been fueled by the infall of fresh metal-poor gas. This would have altered drastically the AMR, with new stars forming at lower metallicity and resulting in a comparatively blue HB morphology.

At present, it is not clear how to identify whether more distant galaxies might exhibit similar behavior, save perhaps from spectroscopic data. However, we note that the incidence in our sample is low, being approximately 10%. Under the assumption that our galaxy sample is broadly representative, which is implicit in the use of our calibration in the first place, this 10% incidence can be considered as a rough estimate of the abundance of such objects.

Eq. 7 allows us to estimate a mass-weighted stellar population age for distant resolved galaxies for which the oldest MSTO is not accessible. The two requirements for using this relation are (i) a CMD that reaches the HB with SNR $>12$ (see § 4.4) and (ii) an estimate of the mean metallicity (e.g., through the LZ relation).

In this section, we validate our results by comparing our HB-based $t^{*}$ values to MSTO-based $t^{*}$ values for Local Group galaxies that are not in our calibration sample. Specifically, we use photometric catalogs obtained from deep HST imaging of the Cetus and Tucana dSph galaxies (HST-GO-10505, PI: Gallart) and of the Eridanus II dwarf (HST-GO-14224, PI: Gallart; HST-GO-14234, PI: Simon; cfr. § 3.1.2). The catalogs have been obtained following the same photometric reduction of our primary sample (Savino et al., 2023). For a consistent comparison, we re-derived the MSTO-based SFHs using the same methodology used for our calibration sample (Savino et al., 2023, Savino et al., in prep.). This choice minimizes systematics due to, e.g., stellar models and and ensures a consistent test of our HB-based calibration. The MSTO-based SFH of Eridanus II is derived from the F606/F814W dataset.

For consistency with the calibration sample, and with application on distant galaxies, we adopted a mean metallicity value from the LZ relation of eq. 2 (Kirby et al., 2013b). For Cetus and Tucana, we used absolute luminosities from McConnachie (2012), while for Eridanus II we used the measurement from Crnojević et al. (2016). Within the respective uncertanties, the metallicities obtained from eq. 2 are in good agreement with more direct metallicity determinations (e.g., Li et al., 2017; Taibi et al., 2018, 2020; Fu et al., 2022). The $M_{V}$ values and the derived $\rm\langle[Fe/H]\rangle$ are reported in Tab. 3.

Figure 7 shows the CMDs of our validation sample and the BHB/RHB selection used to calculate $\eta$. The values of $t^{*}$ inferred from the HB, and those calculated from the MSTO-based SFH are reported in Tab. 3. The agreement between HB-based and MSTO-ages is excellent. In all cases, the two measurements differ by less than 0.4 Gyr, a difference that is small compared to our uncertainty budget ($\sim 0.8$ Gyr). It is also comparable to other common sources of uncertainties in the ages of resolved galaxies (e.g., stellar models, distance and reddening uncertainties, Gallart et al., 2005). The case of Eridanus II also demonstrates consistency between the HB-based ages from F475W/F814W and F606W/F814W data. Despite differences in the exact BHB/RHB selection, highlighted in § 3.1.2, the ages obtained from the two datasets differ by only 0.1 Gyr.

All galaxies analyzed so far in this paper have not experienced significant star-formation activity in the last few Gyr. This means that the number of sources that could contaminate our BHB selection is low. In galaxies with more recent star formation, however, massive stars on the MS will overlap in magnitude and color with the stars on the blue end of the HB, complicating the measurement of the HB morphology. In this section we apply our methodology to two such galaxies, to quantify the impact of this effect.

Figure 8 shows the F475W/F814W CMD of the Pisces and Pegasus dwarf irregular galaxies, which are part of the M31 satellite Treasury sample. The CMDs clearly show the presence of a plume of bright MS stars in our BHB selection box. We deal with this contamination using the same method described in § 3.1.3. We determine the luminosity function of the MS using stellar counts above and below the BHB box (magenta boxes in Fig. 8). The MS boxes span 0.5 mag in magnitude and are subdivided in five and three bins for the lower and upper boxes, respectively. The number of MS contaminants, $n_{MS}$, is then estimated by fitting the luminosity function with a power-law and integrating the resulting luminosity density over the span of the BHB box. The decontaminated sample of BHB stars is then used to calculate $\eta$.

The values of $t^{*}$ we measure through the HB morphology of Pisces and Pegasus are $11.48^{+0.70}_{-0.66}$  Gyr and $10.31^{+0.74}_{-0.69}$ Gyr, respectively (assuming absolute luminosities from Savino et al. 2022). We compare these values with the corresponding $t^{*}$ obtained from the MSTO (Savino et al., in prep., listed in Tab. 3) and find that the two measurements differ by approximately 0.6 Gyr. This is somewhat larger than the difference measured for passive galaxies in § 4.1, which reflects the added uncertainty introduced by the MS decontamination. Nevertheless, HB and MSTO based ages are still compatible to the level of both statistical and systematic uncertainties. The conclusion is that, while the presence of the bright MS might have an effect on the fidelity of HB-based ages, the accuracy of our method is still marginally impacted.

The ability to obtain reliable stellar ages from the morphology of the HB is of particular value for galaxies in the immediate vicinity ($1\lesssim D\lesssim 5$ Mpc) of the Local Group. Detecting the oMSTO of galaxies in this distance range is either impossible or impractically expensive with current telescopes. On the other hand, the brighter and less crowded HB stars can be observed much more efficiently in these galaxies and represent the best available tracer of their early star formation. In this section, we provide a first application of our new method by deriving stellar ages for the isolated galaxies VV 124 ($D\sim 1.37\pm 0.07$ Mpc, Tully et al., 2013) and KKR 25 ($D\sim 1.93\pm 0.06$ Mpc, Makarov et al., 2012).

VV 124 and KKR 25 are two relatively bright ($10^{6}\lesssim L/L_{\odot}\lesssim 10^{7}$), metal-poor ($-1.9\lesssim\langle[Fe/H]\rangle\lesssim-1.5$) dwarf galaxies located just outside the LG (e.g., Karachentsev et al., 2001; Kopylov et al., 2008; Bellazzini et al., 2011b; Makarov et al., 2012; Kirby et al., 2012, 2013a). Both galaxies host predominantly old/intermediate stellar populations with a small amount of recent $(\lesssim 1$Gyr) star formation (Kopylov et al., 2008; Bellazzini et al., 2011a; Makarov et al., 2012; Neeley et al., 2021), have relatively low $\rm H_{I}$ content (Begum & Chengalur, 2005; Bellazzini et al., 2011b), and (in the case of VV 124) pressure-supported kinematics (Kirby et al., 2012). Many of these observations have led to suggestions that these galaxies are transition dwarfs in a late stage of their evolution.

The properties of KKR 25 and VV 124 are particularly intriguing in the context of low-mass galaxy evolution. The long-established morphology-density relationship has been central to the idea that environment is a primary driver in the evolution and quenching of luminous ($L\gtrsim 10^{6}L_{\odot}$) dwarf galaxies (e.g., Mateo, 1998; Grebel et al., 2003; Geha et al., 2012; Putman et al., 2021). Yet, early-type dwarfs such as KKR 25 and VV 124, are distant enough from other massive galaxies that past interactions are virtually impossible. This poses a challenge in our understanding of why such objects have relatively low levels of recent star formation and small neutral gas reservoirs in light of environmentally-driven theories of galaxy quenching.

While the distances of these two galaxies make their oMSTOs very challenging to observe, their HBs are well within the reach of existing HST observations. To measure the HB morphology of these two targets we have therefore performed PSF-photometry (using the same procedure as in § 2.1) on the archival HST data of program GO-15244 (PI: Monelli). Each galaxy has 16 orbits of ACS F475W and F814W imaging that have previously been used to study their RR Lyrae (Neeley et al., 2021).

Figure 9 shows the CMDs of KKR 25 and VV 124 around the HB region. The helium burning sequences of these galaxies are extended from the blue end of the HB to more massive helium-burning stars above the red clump, suggestive of an extended SFH. A notable feature in this figure is that both galaxies have a bright plume of young main-sequence (MS) stars which overlaps with the BHB. We account for this using the same procedure described in § 4.2.

The resulting HB morphology parameters are $\eta=0.063\pm 0.014$ for KKR 25 and $\eta=0.132\pm 0.007$ for VV 124. Using absolute magnitudes from McConnachie (2012) and deriving average metallicities through eq. 2 (listed in Tab. 3), we find $t^{*}=11.21^{+0.69}_{-0.65}$ Gyr for KKR 25 and $t^{*}=11.03^{+0.73}_{-0.68}$ Gyr for VV 124. These mass-weighted mean ages are similar to Cetus and Pisces, and quite younger than Tucana and Eri II.

Given that $t^{*}$ is an integrated measurement of the SFH, the interpretation of the mean ages of KKR 25 and VV 124 depends on our assumptions on when star formation began in these objects. Under the assumption that star formation began at the earliest epochs, our age measurements indicate a substantially extended SFH, akin to the aforementioned Pisces and Cetus dwarfs (Monelli et al., 2010; Hidalgo et al., 2011). On the other hand, it is possible that star formation remained low for the first one or two Gyr, and then experienced a short duration burst around $z\sim 3$. While we cannot discard this option, we note that VV 124 presents a prominent BHB population (clearly seen through the MS contamination), which suggest the presence of a sizable population of truly ancient stars.

Regardless, the presence of significant star formation at $\sim 11$ Gyr suggests that reionization had little impact on the early star formation of these galaxies. In comparison, isolated galaxies Leo A and DDO 210 have similar present day luminosities, but only formed 10% of their stellar mass at ages $>$10 Gyr, which, in some interpretations, could be due to suppression from the cosmic ultra-violet background (e.g., Cole et al., 2007, 2014). Similarly, the degree of isolation of these galaxies, combined with their star formation activity at low ($\lesssim 3$) redshift, suggests that the transition from gas-rich to gas-poor dwarf galaxies cannot solely be a function of environment. According to the 2023 version of the Nearby Galaxy Catalog (Karachentsev et al., 2004, 2013), the nearest massive neighbors to KKR 25 and VV 124 are more than 1 Mpc away (UGC8508 and Sextans B, for KKR 25 and VV 124, respectively), and even those are relatively low-mass dwarfs ($M_{*}\sim 10^{8}M_{\odot}$). Any environmental interaction these galaxies might have had with a large galaxy must have happened at ancient times and, if responsible for their nature of early type galaxies, would have resulted in a very early depression of star formation.

For reference, the Cetus dSph, which has been proposed as a candidate backsplash galaxy from the MW or M31 (e.g., Sales et al., 2007; Teyssier et al., 2012), has a similar $t^{*}$ to VV 124 and KKR 25 but it currently lies at $\sim 700$ kpc from either large spirals (McConnachie, 2012). This is roughly two times closer than VV 124 and three times closer than KKR 25. Tucana, another candidate backsplash galaxy, lies at a larger distance of $\sim 900$ kpc from the MW (McConnachie, 2012), but has a $t^{*}$ that is $\sim 1$ Gyr older. It therefore seems unlikely that VV 124 and KKR 25 could have had an interaction with the MW or M31 and managed to sustain star formation for as long as their HB seems to indicate.

Finally, if, these galaxies had ancient SFHs that are more like MW satellites, i.e., dominant ancient episodes, than other isolated galaxies (e.g., Leo A, DDO 210), then they provide interesting counter examples to the slow/fast scenario proposed by Gallart et al. (2015), in which the dominance of ancient star formation can be mapped to the local environmental density at high-redshifts. While this limited sample hints at new insight into isolated dwarf galaxy formation, a large systematic study (e.g., if HB-depth photometry were available for the entire ANGST sample; Dalcanton et al. 2009) is needed to enable a truly statistical study.

A challenge of characterizing the HB morphology of galaxies beyond the LG will be the limited photometric precision that observations of these distant targets will achieve at the HB magnitudes. In contrast, our calibration sample has exquisite photometry of the HB region, with mean SNR ratios between 147 (And XVI) and 58 (And XXVI). To quantify the effect of lower photometric precision on the recovered ages, we have artificially degraded the photometry of one of our calibration galaxies and re-applied our methodology.

Figure 10 illustrates our experiment. We start with the original photometry of And XIV as our ground truth. The mean magnitude SNR of this dataset, at the magnitude level of the HB, is 74. The measured HB morphology index is $\eta=0.29\pm 0.04$, which results in an age of $12.12\pm 0.71$ Gyr. We then perturbed the photometry of individual stars in And XIV with a random Gaussian term $\mathcal{N}(0,\sigma)$. The term $\sigma$ is chosen to emulate a target SNR: $\sigma=\sqrt{\frac{1}{SNR^{2}}-\frac{1}{74^{2}}}$. We adopt the same $\sigma$ for the whole photometry, as we are considering a relatively small magnitude range in the CMD. The upper panel of Fig. 10 illustrates select examples of the CMDs we obtain. We then measure $\eta$ from the new degraded CMDs and derive $t^{*}$ values using eq. 7.

Because the morphology of the HB is parameterized using large selection boxes, the value of $\eta$ only begins to significantly change as a function of SNR for SNRs as low as $\sim 10$. We find that below SNR$\sim 10$, the measured value of $\eta$ can deviate significantly from the truth, in excess of measurement uncertainties. Though the $t^{*}$ is not heavily affected in our example, because the scatter term in Eq. 8 dominates, we find that the shallowest photometric depths affect our accuracy to $\sim 200$ Myr. Moreover, at such low SNRs, the lower RGB will no longer be accessible, meaning that our decontamination will not work. In contrast, for HBs with SNR $\gtrsim 10$, $\eta$ and age are statistically consistent and the lower RGB can be accessed for decontamination purposes.

Thanks to future deep, large-area surveys (e.g., Rubin, Roman, Euclid), our census of low-mass, nearby galaxies is poised to dramatically increase in the coming decade. For example, a wealth of new low-mass galaxies down to $M_{V}\sim-5$ within $\sim 5$ Mpc should be discoverable (e.g., Mutlu-Pakdil et al., 2021). Characterizing the ages of these galaxies from the MSTO will be challenging due to the large distances. While the HB will be much brighter, it may also be sparsely populated (e.g., a few dozen HB stars for bright ultra-faint dwarfs). Accordingly, it is important to quantify the accuracy of our method as a function of stellar population size.

To do so, we start by using two of our calibration galaxies, with different values of $\eta$, for ground truth. We have selected And V ($M_{V}=-9.3$; Savino et al. 2022), which has a population of $n_{red}+n_{blue}=1651$ stars (before decontamination) and a relatively red HB morphology ($\eta=0.14$), and And XVII ($M_{V}=-7.8$; Savino et al. 2022), which has a population of $n_{red}+n_{blue}=479$ stars and a much bluer HB morphology ($\eta=0.72$). We then re-sample (without replacement) the CMD of these two galaxies until a given number of HB stars, $n_{HB}=n_{red}+n_{blue}$, is reached. For a given $n_{HB}$ we generate 1000 CMD realizations and calculate the corresponding $\eta$ and $t^{*}$ values. We also scale the absolute luminosity of both galaxies accordingly, as we decrease the HB counts, to estimate the approximate galaxy luminosity at which HB sampling becomes relevant.

Figure 11 shows the difference between the original measurements of $\eta$ and $t^{*}$, and those measured from the re-sampled CMDs, as a function of $n_{HB}$. We do not report any significant bias in the measured ages, even in sparsely populated CMDs. This is likely because stochastic sampling affects stars equally regardless of their position on the HB. The reduced HB population, however, does steadily increase the measurement uncertainty on $\eta$, as it is reasonably expected. Nevertheless, the bottom panels of Fig. 11 show that the contribution of stochastic sampling to the uncertainty on $t^{*}$ remains subdominant with respect to our calibration scatter for HB populations of $\gtrsim 70$ stars for And V and $\gtrsim 30$ stars for And XVII. Below these HB population regimes, the low number of HB stars becomes the primary source of uncertainty in age.

The effect of low HB population kicks in at different HB counts for And V (70 stars) and And XVII (30 stars) due to their different HB morphology. Given that eq. 7 depends on the logarithm of $\eta$, galaxies with blue HB morphologies (i.e., $\eta\sim 1$) will be affected less by a given change in $\eta$ than galaxies with red HBs (i.e., $\eta\sim 0$). The corresponding absolute luminosity at which HB sampling becomes relevant is $M_{V}\sim-6$ for And V and $M_{V}\sim-5$ for And XVII. The galaxy luminosity limit for our method applicability is therefore slightly dependent on HB morphology. We note, however, that at such low magnitudes, blue HBs are typically more common.

The above experiment assumes that the only source of uncertainty introduced by a low HB population is due to stochastic sampling. In reality, at low HB counts the issue of identifying bona fide HB stars become increasingly more relevant. Spurious sources in the HB region, such as photometric artifacts, foreground stars, or background galaxies, can potentially bias the inferred value of $\eta$ and, consequently, of $t^{*}$. Such an effect is much more difficult to quantify, as it depends on a number of other factors, such as the depth of the photometry, the distance of the galaxy, its galactic latitude, its SFH, and the adopted methodology to reject non-stellar sources. As a general rule, we advise using our calibration with caution whenever the number of potential contaminants, as estimated by visual inspection of the CMD, is a non-negligible fraction of the stars in the HB region.