A detailed understanding of stellar evolution remains one of the most sought-after goals in astrophysics. Popular stellar evolution codes such as Geneva (Schaller et al. 1992;Schaerer et al. 1993), Y 2 (Kim et al. 2002;Yi et al. 2003;Demarque et al. 2004), BaSTI (Pietrinferni et al. 2004(Pietrinferni et al. , 2006;;Hidalgo et al. 2018), Dartmouth (Dotter et al. 2008), PARSEC (Bressan et al. 2012;Chen et al. 2014), and MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018) ) are powerful tools for interpreting observations of stellar populations. However, many essential topics in stellar evolution are still not well understood and/or not properly taken into account in even the most state-of-the-art models. Examples of such topics include convection, mass loss and mass transfer, common envelope evolution, and binary interaction. Often, these three-dimensional phenomena are approximated by simplified parametric models tuned to specific observables and integrated into one-dimensional stellar evolution codes. These uncertainties limit our understanding of many important phases of stellar evolution, such as the horizontal branch, the asymptotic giant branch (AGB) and post-AGB phase, planetary nebulae, and supernovae (see discussions in Gallart et al. 2005;Conroy et al. 2009;Conroy 2013).
The delay-time distribution (DTD) is a promising method for testing stellar evolution models in complex stellar populations (Maoz & Mannucci 2012;Maoz et al. 2014). The DTD is defined as the occurrence rate of a class of astronomical object as a function of time since a hypothetical brief burst of star formation; it is equivalent to the impulse response, or Green's function. The DTD constrains the evolutionary timescale and formation efficiency of the object's progenitors, and theoretical DTDs are common predictions of stellar population synthesis models (Mennekens et al. 2010;Nelemans et al. 2013;Toonen et al. 2013;Zapartas et al. 2017). Observationally, DTDs can be derived from surveys of objects, provided that the stellar age distributions (SADs) of their host galaxies are measured (Gal-Yam & Maoz 2004;Totani et al. 2008;Maoz & Sharon 2010;Maoz et al. 2011Maoz et al. , 2012;;Graur et al. 2014;Maoz & Graur 2017;Friedmann & Maoz 2018).
More recently it was shown that DTDs can be a useful diagnostic of stellar evolution in Local Group galaxies with highquality observations of their resolved stellar populations (Badenes et al. 2010(Badenes et al. , 2015;;Maoz & Badenes 2010). Using this approach, Badenes et al. (2015) measured the first DTD for planetary nebula, showing evidence of two distinct populations of planetary nebula progenitors: one with ages of 5-8 Gyr and another with ages of 35-800 Myr. The key advantage of a DTD is that it constrains the evolutionary timescales for the progenitors of the entire population of objects in a galaxy, taking into account the variety of star formation histories that these objects have evolved from. It can be a powerful tool for identifying or ruling out the presence of specific formation channels, measuring their efficiency, and identifying physical mechanisms that are not part of canonical progenitor models.
In this paper, we will use a SAD map of the Large Magellanic Cloud (LMC) to derive the DTD of RR Lyrae stars-pulsating horizontal branch stars with periods between 0.2 and 1 day (Smith 2004). We chose to test the DTD method on RR Lyrae for several reasons. First, the sample size of RR Lyrae in the LMC is quite large (see Section 2.1), allowing us to measure a DTD with high significance. Second, there is strong evidence that RR Lyrae are mostly ancient stars, older than 10 Gyr, given their pulsational properties (Smith 2004;Marconi et al. 2015) and abundance in old globular clusters (Clement et al. 2001;Soszyński et al. 2014Soszyński et al. , 2016)). Measuring an RR Lyrae DTD therefore provides a rigorous test of the DTD method for the recovery of progenitor age distributions as well as measurements of the star formation history of old resolved stellar populations. Lastly, a DTD analysis provides an opportunity to test stellar evolution models of RR Lyrae in a completely new way. While there is consensus on the interpretation of RR Lyrae as ancient stars, the lower limit on their ages has been somewhat unconstrained. For example, the absence of RR Lyrae in the SMC cluster Lindsey 1 ( » t 9 Gyr), compared to their presence in NGC 121 ( » t 10-11 Gyr), is generally cited as evidence of a lower limit of 10 Gyr for the progenitor age of RR Lyrae (Olszewski et al. 1996;Glatt et al. 2008). However, growing evidence of thin-disk RR Lyrae in the Milky Way raises questions about whether this limit might be lower, and whether an intermediate-age progenitor channel can exist (Layden 1995;Zinn et al. 2019;Prudil et al. 2020). Additionally uncertainties may exist in the evolutionary models of RR Lyrae stars themselves: since RR Lyrae are horizontal branch stars, their positions on the color-magnitude diagram depend on an unknown combination of factors (commonly known as the second-parameter phenomenon) such as metallicity, age, mass loss on the red giant branch, stellar rotation, core structure, and chemical abundance (see Fusi Pecci & Bellazzini 1997;Catelan 2009;Dotter 2013, for reviews).
The paper is organized as follows. Section 2 briefly describes the two ingredients for calculating the RR Lyrae DTD-the OGLE-IV survey of RR Lyrae, and the SAD map of the LMC from Harris & Zaritsky (2009). Section 3 describes the measurement of the RR Lyrae DTD from the OGLE-IV survey. Section 4 checks whether the measured DTD is consistent with the observed spatial distribution of RR Lyrae in the OGLE-IV survey, incompleteness in the SAD map, and RR Lyrae in star clusters. Section 5 discusses the two possible interpretations of the DTD result-that RR Lyrae may have a previously unidentified younger progenitor channel, and that unknown systematics still exist in SAD measurements of resolved stellar populations older than a gigayear.
We use the catalog by Soszyński et al. (2016) of 39,082 RR Lyrae stars from the Optical Gravitational Microlensing Experiment (OGLE-IV)foot_0 survey of the LMC (Udalski et al. 2015). These RR Lyrae were selected by the OGLE-IV pipeline from the full database of OGLE I-band light curves with periods between 0.2 and 1 day, and then further classified as fundamental (RRab), first-overtone (RRc), and mixed-mode pulsators based on their periods, Fourier amplitudes, and lightcurve shapes. We excluded catalog entries that were flagged as Galactic RR Lyrae or eclipsing variables, objects with uncertain classification, and sources that fall outside the SAD maps. This produced a final sample of 29,810 RR Lyrae (Figure 1). The photometric completeness of the sample is quite high, as evidenced by the I-band luminosity function of RR Lyrae in the most crowded OGLE-IV field (LMC 503, see Figure 2). The field has a completeness limit of » I 20.5 (Udalski et al. 2015), whereas the RR Lyrae sample in the field has a median magnitude ¯= I 18.82 and standard deviation s = 0.4 I . The median magnitude of the RR Lyrae is nearly s 4.2 I above the completeness limit (almost 99.9% of the RR Lyrae have I-band magnitude above the completeness limit). Thus even in the most crowded region, the RR Lyrae sample can be considered photometrically complete. Also as seen in Figure 2, the I-band luminosity function of the full OGLE-IV sample inside the region of Harris & Zaritsky (2009) (see Section 2.2) is also well above the completeness limit of the most crowded region. Such high completeness is important for measuring unbiased rates and DTDs (Badenes et al. 2010;Maoz & Badenes 2010).
We use the SAD map of the LMC constructed by Harris & Zaritsky (2009, hereafter HZ09). This map provides the best-fit stellar mass formed as a function of lookback time in spatial cells resolving the central ´ 8 . 5 7 . 5 of the galaxy (Zaritsky et al. 2004). HZ09 also provides the associated 1σ upper and lower limits to the SAD in each cell, which we will incorporate into the DTD uncertainties in Section 3.2. The SADs were calculated using data from the Magellanic Cloud Photometric Survey (MCPS) of nearly 4 million stars collected with the 1 m Swope telescope, down to a completeness of V=20-21 mag (Zaritsky et al. 1997(Zaritsky et al. , 2004)). The MCPS region was divided into 1376 cells, each measuring ¢ ´¢ 24 24 , or ¢ ´¢ 12 12 if the field contained more than 25,000 stars. The contours of these cells are shown in Figure 1. SADs were derived using the StarFISH algorithm by fitting the color-magnitude diagram in each cell with a linear combination of isochrones (Harris & Zaritsky 2001). After accounting for extinction and photometric errors, each cell's SAD was fit using 16 logarithmically spaced bins spanning the ages between 4 Myr and 20 Gyr, and four metallicity bins (Z=0.008, 0.004, 0.0025, and 0.001). For ages younger than 100 Myr, a single metallicity of Z=0.008 was used because the different metallicity isochrones are almost indistinguishable. Note that while the lower limit on the age of RR Lyrae stars is generally quoted as 10 Gyr, the SAD map has a single indivisible age bin of 8-12 Gyr, and so we refer to this lower limit as 8 Gyr in the rest of the paper.
The RR Lyrae catalog and SAD maps of the LMC are used to estimate the RR Lyrae DTD using the non-parametric method described in Badenes et al. (2015, hereafter B15), although we improve on some aspects of it. The RR Lyrae DTD is the number of RR Lyrae formed per unit stellar mass as a function of the time delay between star formation and the RR Lyrae phase. The convolution of the DTD with the SAD in each spatial cell predicts the number of RR Lyrae that will be produced by the stellar population in that cell. For each SAD cell i, the number of RR Lyrae predicted (λ i ) is
where M ij is the stellar mass formed in age bin j and cell i, Ψ is the RR Lyrae formation rate (RR Lyrae per year per unit M e ), and T vis is the duration of the RR Lyrae phase (both Ψ and T vis are functions of age bin). The widths of the age bins, j, are selected using the methodology described in Appendix A and provide the best compromise between detection significance and temporal resolution.
In this paper, we retain the notation of B15 and refer to the quantity ( ) YT vis as the DTD (with units of RR Lyrae per unit M e ). The DTD (ΨT vis ) is determined by minimizing the difference between the predicted and observed numbers of RR Lyrae across spatial cells. This is carried out in a similar manner to B15 using the Markov Chain Monte Carlo (MCMC) solver emcee (Foreman- Mackey et al. 2013). We denote [ ] = N N i as the vector representing the number of RR Lyrae in each spatial cell, and
as the vector of the predicted number of RR Lyrae per stellar mass for each age bin . The posterior is calculated as
where is the likelihood and π is the prior. We assume is a product of either Poisson or Gaussian probabilities depending on N i ,
where K is the number of SAD cells. Although Poisson distributions converge to Gaussian for large N i , we specify the likelihoods separately to avoid computational issues. The prior is defined such that the ( ) YT j vis values have uniform probabilities in log-space:
We propagate the errors on the M ij values in the SAD map into uncertainties on YT vis . B15 calculated DTDs for the bestfit SAD, the upper limit on the SAD, and the lower limit on the SAD. The differences between the best-fit and upper/lower limits were treated as the DTD's 1σ uncertainties.
In this paper, we use an improved method of propagating SAD uncertainties into the DTD. We randomly generate 100 mock SADs, assuming M ij has the normally distributed uncertainties given by HZ09, and calculate a DTD for each mock SAD. We combine the MCMC posterior chains from these 100 DTDs into a single chain, and estimate the 95% credible interval on this chain using a criterion of highest posterior density. We define the mode of the distribution minus the upper and lower limits of this interval as our 2σ + and 2σ - confidence intervals, respectively. We define a "signal" detection in each bin j of the DTD as a value of ΨT vis in that bin that is 2σ -above 0. Non-detections are presented as 2σ limits on the DTD signal in a particular age bin.
We estimate the contribution of each age bin to the total sample by multiplying the value of ΨT vis by the total stellar mass formed in each age bin. The percentage contribution will then have uncertainties due to both the DTD and the total stellar masses formed. We estimate these uncertainties with a Monte Carlo method. We draw DTD values from the recovered posterior probability distributions (called ( ) YT j vis ), multiply by the total stellar mass per age bin in each randomized SAD map (called M j ), and get the number of RR Lyrae contributed by age bin, l j . The percentage contributions ( / l l = å j j ) and their 1σ uncertainties are listed in Table 1.
The DTD for the full OGLE-IV RR Lyrae sample and the two main RR Lyrae subtypes is shown in Figure 3, and the values are tabulated in Table 1. We detect a signal in the DTD at a high significance (>5σ) for all age bins older than 1.3 Gyr. The detections with the highest significance are found in the 2.0-3.2, 5.0-7.9 and >12.6 Gyr bins. Most of this signal is contributed by the RRab stars, which are the most common subtype in the OGLE-IV sample. The RRc subclass contributes to the DTD mostly above 2 Gyr. Although RRc's have been susceptible to confusion with other variable sources in timedomain surveys (Kinman & Brown 2010;Mateu et al. 2012;Drake et al. 2014), Figure 3 shows that the total DTD is dominated by RRab objects in all age bins, making it unlikely that the DTD is biased by sample contamination. The detected signal of the full DTD is relatively flat above 1.2 Gyr, with about 1-3 RR Lyrae produced per 10 5 M e of stellar mass formed. While there is a signal in the DTD below 0.8 Gyr, it falls below our 2σ -threshold, and we only show the 2σ + upper limit.
Our detection of a DTD signal below 8 Gyr is a surprising result. About 46% of the LMC's RR Lyrae stars are produced from populations older than 10 Gyr, the age range generally inferred for RR Lyrae in star clusters. But about 51% of the DTD signal comes from progenitors with ages between 1.3 and 8 Gyr, and this result has high (>5σ) significance. If we recalculate the DTD without assuming normal errors on M ij (i.e., using the same method as B15, see Appendix B), we still detect a strong signal below 8 Gyr, but with a total contribution of 41%. A comparison of these timescales with those of the PARSECfoot_2 models for the onset of helium burning leads to the conclusion that RR Lyrae can arise from main-sequence progenitors as massive as ∼2 M e at LMC metallicity. Incidentally, this upper limit is similar to the mass at which stars transition from igniting He under degenerate conditions (the "He flash") to burning helium under stable, non-degenerate conditions (Bildsten et al. 2012;Mosser et al. 2014). We also tested the dependence of the DTD on the period and brightness of RR Lyrae stars by subdividing the OGLE-IV data by pulsation time (three bins with <0.45day, 0.45-0.58day, and >0.58day) and I-band magnitude (three bins with I<19. 2, 19.2-19.4, and >19.4). By subdividing the sample in this way, we ensured that each period and magnitude bin contained enough RR Lyrae for a robust measurement of their DTDs. We found that all the subsamples have significant (>2σ) DTD signals in the range 1-8Gyr, with no discernible trend in the shape of the DTD with magnitude or period.
Although, in principle, it is possible to derive a metallicitydependent DTD with the HZ09 SAD maps, we defer that study to a future work. However, it is well known that LMC RR Lyrae are generally metal-poor, with [Fe/H]<-0.5 and a peak in the metallicity distribution function at [Fe/H]∼-1.5 (Skowron et al. 2016). We can therefore check to see if there is DTD signal below 8 Gyr under the assumption that OGLE-IV RR Lyrae are only produced by the HZ09 SAD in the two metal-poor bins, i.e., for < Z 0.0025 or [Fe/H]-1.02 (Bertelli et al. 1994). Star formation in this metallicity range dominated the LMC SAD until ∼2 Gyr, as seen in Figure 4(a), and the resulting DTD is shown in Figure 4(b). Even if we restrict our analysis to the metal-poor SAD, the signal below 8 Gyr persists. This implies that, even if we assume that all LMC RR Lyrae are produced by metal-poor stars, progenitors younger than 8 Gyr are still needed to explain their distribution.
Before delving into the implications for RR Lyrae studies, we carry out a more detailed analysis of the recovered DTD to assess the robustness of our results.
The HZ09 SAD map of the LMC allows us to test whether the traditional DTD for RR Lyrae (i.e., one in which the progenitors are always older than 8 Gyr) is consistent with the spatial distribution of RR Lyrae observed by OGLE. This is equivalent to inverting the DTD recovery process. To do this, we generate mock RR Lyrae maps by convolving the SAD map with a DTD that is non-zero only in the two oldest age bins, 8-12 and 12-20 Gyr. The total stellar mass that formed in these age bins was 1.23 10 9 M e . Assuming all 29,810 RR Lyrae were produced by progenitors in these age bins results in a DTD of the form We generate 100 mock RR Lyrae maps using this DTD, where the number of RR Lyrae per cell, N i , is drawn from the Poisson distribution in Equation (3) with λ i given by Equation (1). Figure 5 shows the DTD recovered from this analysis. Our MCMC solver correctly measures the input mock DTD with strong detections only in the two oldest age bins. The younger bins have 2σ upper limits that are almost an order of magnitude lower than the DTD recovered from the OGLE-IV sample. Moreover, the difference between populations of RR Lyrae with the OGLE-IV and mock DTDs can be seen visually in Figure 6. The measured OGLE-IV DTD predicts a distribution of RR Lyrae stars elongated along the LMC Bar and declining smoothly with radius (center bottom panel of Figure 6); this is very similar to what is seen in the OGLE-IV map (left panel). In contrast, our mock distribution from a uniformly old RR Lyrae progenitor population (Equation ( 5)) has more structure and leaves larger residuals in the difference map (top center and right panels of Figure 6). To identify significant regions of discrepancy in the OGLE-IV DTD and mock old DTD maps, we generate 10 4 maps of RR Lyrae, where in each map the number of RR Lyrae in each cell is drawn from a Poisson distribution with mean number of RR Lyrae in that cell given by the DTD according to Equation (1). We then estimate the mean and standard deviation of the number of RR Lyrae per cell, and identify cells with black squares in Figure 6 where the observed number of RR Lyrae is greater than five times the standard deviation from the mean. We find the mock old DTD map has a larger number of cells where the RR Lyrae number count is discrepant, and these cells are mostly located in the Bar region, with a few located outside (note that this is only for the purpose of visualizing the discrepancy; the actual DTD is constrained by the joint likelihood measured from all the cells as in Equation (3)).
We conclude that if the OGLE-IV RR Lyrae had indeed been produced exclusively from the old ( > t 8 Gyr) stars of the HZ09 SAD map, our method would have recovered the correct DTD, without spurious signals at younger ages. The DTD signal we recover at ages between 0.8 and 8 Gyr must either be real or be an artifact produced by systematics in the HZ09 SADs. Next, we investigate the impact of these SAD systematics on the DTD more thoroughly.
We explored whether any incompleteness in the MCPS photometry (from which the SAD was measured) could be driving the intermediate-age signal. The photometric completeness limit is around V=21 mag, with the Bar region being almost 1 mag shallower as a result of the stellar crowding. Because of this, StarFISH could reliably solve for the number of stars per age bin only for ages younger than 4 Gyr (although we note that our intermediate-age DTD signal is detected younger than 4 Gyr). Populations older than 4 Gyr were traced by their giant branch stars, and assigned by StarFISH to a single age bin. To anchor the "shape" of the SAD beyond 4 Gyr, stellar ages were measured in a few isolated Hubble Space Telescope (HST) fields in the LMC Bar (Holtzman et al. 1999;Olsen 1999;Smecker-Hane et al. 2002). As shown in Figure 7 of HZ09, the SADs of each HST Bar field are characterized by star formation at lookback times of 10 and 5 Gyr, with a period of quiescence in between. This consistency allowed HZ09 to adopt a common shape for the SAD in the rest of the Bar region.
We first checked whether an underestimation of the total mass formed at old stellar ages due to any photometric incompleteness could be driving the intermediate-age signal. Obtaining a complete census of old stellar mass is generally difficult without including infrared data (Conroy 2013), so we studied the changes to our DTD by manually changing the old SADs. We recalculated the RR Lyrae DTD with the same SAD map, except with each cell's 8-12 and 12-20 Gyr stellar mass multiplied by factors of 2, 4, and 10. We find that regardless of how much mass is added, the DTDs still find signals for ages between 0.8 and 8 Gyr with >5σ confidence. Since the incompleteness in the SAD map is dominant in the Bar region, we also tried recalculating the DTD without the LMC Bar (i.e., by removing the SAD cells and their RR Lyrae from our calculations). We show the RR Lyrae DTDs without the "Inner Bar" in Figure 7, and without both the "Inner Bar" and "Outer Bar" in Figure 8 (the cells of these regions are defined according to Figure 6 in HZ09). The recovered DTDs outside the excluded regions are similar to the DTD obtained from the full OGLE-IV sample. For the case with both the Inner and Outer Bars removed, the number of RR Lyrae and the number of SAD cells are smaller, leading to a higher upper limit on the DTD older than 0.8 Gyr, and a non-detection in the 3-5 Gyr bin. However, there is still significant signal at ages younger than 8 Gyr, and particularly below 4 Gyr, where the ages are determined by the main-sequence turnoff in the MCPS photometry. It therefore appears unlikely that missing stellar masses at ages >8 Gyr due to incompleteness or to crowding in the Bar region are the only factors driving the intermediate-age signal in the DTD.
Our measured DTD suggests that progenitors of RR Lyrae can be as young as ∼1 Gyr, which is much smaller than the lower limit on the RR Lyrae age of 10 Gyr inferred from star clusters (Olszewski et al. 1996). We therefore examine the number of OGLE-IV RR Lyrae in the LMC star clusters as a function of age, realizing that cluster membership can only be confirmed with spectroscopic and proper-motion measurements that are beyond the scope of this work. . Bica et al. (2008) note that these are "apparent" sizes (measured as far as the background limit in the images), but they most likely enclose the majority of the cluster mass. We also obtained chemical abundances for 29 clusters ([Fe/H], shown in Table 2 with references), and as we show later, this is sufficient for understanding the correlation of cluster RR Lyrae statistics with metallicity. Although there have been more recent studies of the population of LMC star clusters (e.g., Nayak et al. 2016;Bitsakis et al. 2017;Piatti 2017Piatti , 2018)), the effects of field-star contamination and asterisms on cluster identification in these Out of the 296 clusters, we investigate the RR Lyrae statistics of 62 that have ages above 1 Gyr in order to compare with the DTD. We define n o as the number of OGLE-IV RR Lyrae observed within a circle of radius r c centered on the cluster, and n b as the expected number of RR Lyrae within the cluster area but unassociated with the cluster (we refer to these as "background" RR Lyrae). Assuming N is the number of RR Lyrae within radii r 2 c and r 4 c of the cluster (we allow a buffer region of r c to account for uncertainty in the actual extent of the cluster), we define n b as
Values for n o and n b are listed in Table 2 and shown in Figure 9. LMC globular clusters older than 10 Gyr clearly host populations of multiple RR Lyrae that exceed the background. The richest population is in NGC 1835 with over 100 RR Lyrae, while NGC 2005, NGC 1928, and NGC 1939 p , the difference between the background-subtracted number of RR Lyrae observed per cluster and the number predicted by the DTD given the cluster mass and age (see Equations ( 6), ( 7), and Section 4.3 for details). Colors indicate the cluster mass. The inset zooms in on clusters between 1 and 4 Gyr. Right: Δn vs. metallicity for a subset of 29 clusters with metallicity information available in the literature. The color bar indicates the cluster age. The inset zooms in on the 1-4 Gyr clusters, which also happen to cluster in the metallicity range -1[Fe/H]-0.25. also the case in the Milky Way, where studies of variable stars in well-known intermediate-age systems such as M67 (Pribulla et al. 2008), NGC 188 (Zhang et al. 2004), NGC 6791 (de Marchi et al. 2007), and NGC 6253 (Kaluzny et al. 2014) have found no RR Lyrae that are likely to be cluster members. In contrast, the Milky Way's globular cluster system contains almost 2000 such objects (Clement et al. 2001;Soszyński et al. 2014).
While this deficit of RR Lyrae in intermediate-age clusters appears contradictory to the 1-8 Gyr signal in the DTD, we note that the intermediate-age clusters in the LMC are generally less massive than old globular clusters, and are therefore less likely to host relatively short-lived objects such as RR Lyrae stars. We quantify this with n p , the expected number of RR Lyrae per cluster given our measured DTD ΨT vis , the mass M c , and age t c of a cluster:
For each cluster, we measure
p , i.e., the difference between the observed background-subtracted number of RR Lyrae (n on b ) and the number predicted by the DTD. We estimate the uncertainty in Δn in the following way: we create 10 6 random samples of the 62 clusters, each sample having the same n o per cluster as given in Table 2 but with n b and n p generated from a Poisson distribution with mean n b and n p values given in Table 2. For any cluster in these samples with a negative value of n on b , we set that value to zero. The mean and standard deviation of this random sampling are taken to be the value and error of Δn. These Δn are shown as a function of cluster age and metallicity in Figure 9. We see that Δn is consistent with zero within 2σ for most intermediate-age clusters in the LMC. This is because the DTD predicts ∼1-3 RR Lyrae per 10 5 M e of stars, whereas the intermediate-age clusters have an average mass of 4.3 10 4 M e , so the majority will have less than one RR Lyrae star. This caveat was also raised by Olszewski et al. (1996) in their study of the SMC: if the RR Lyrae-rich cluster NGC 121 ( » t 12 Gyr) were scaled down to the magnitude of the 9Gyr system Lindsey 1 (which has no RR Lyrae), it would host less than one RR Lyrae star. The intermediate-age clusters of the Milky Way, including the ones listed above, have the same issue: all have masses less than 10 4 M e (Kaluzny & Udalski 1992;Mermilliod 2000;Grundahl et al. 2008;Piskunov et al. 2008;Elsanhoury et al. 2016;Kruijssen et al. 2019). Thus, the traditional lower limit on the age of the RR Lyrae progenitor population is subject to small-number statistics.
We note that the globular clusters host a diversity of RR Lyrae populations. NGC 1939, NGC 1916, and NGC 2005 are consistent with Δn=0 within their uncertainties, and their production rates (n on b )/M c =(1.8-2.5)×10 -5 RR Lyrae per M e are similar to the measured DTD. The other old clusters, however, host RR Lyrae significantly in excess of their DTDpredicted n p (i.e., Dn 0), and have an average production rate ∼10 -4 RR Lyrae per M e , a factor ∼5 times higher than the DTD. The high production rate in the old clusters may be an effect of their low metallicities, which is consistent with the differences in production rates of RR Lyrae in the halo versus the bulge and disk of our Galaxy (Layden 1995;Dékány et al. 2018), although the large differences in Δn for the globular clusters could also be an effect of the second-parameter phenomenon in clusters (see Fusi Pecci & Bellazzini 1997;Catelan 2009;Dotter 2013, for reviews). A more detailed assessment of the field and cluster DTDs can be done with a metallicity-dependent DTD that we reserve for future work.
If we assume that the DTD measured from the OGLE-IV survey is correct, there are two possible interpretations of this result: (1) RR Lyrae can form from progenitors younger than 10 Gyr in addition to the conventional route via older stars, and this result was undetected in previous studies due to various observational limitations, or (2) all OGLE-IV RR Lyrae are older than 10 Gyr, and our result is a product of significant (though not readily obvious) systematics in the age derivation of older stellar populations. We discuss both of these interpretations below.
Since the DTD is generally regarded as a reflection of the progenitor age distribution (Maoz & Badenes 2010;Badenes et al. 2015), it is tempting to consider that RR Lyrae in the LMC are being produced by an as-yet undiscovered intermediate-age progenitor channel between 1 and 8 Gyr, in addition to the usual channel older than 8 Gyr.
However, looking at currently available evidence, the possibility of such an undiscovered intermediate-age channel appears to be questionable. Although the small-number statistics described in Section 4.3 is a factor, it is nevertheless true that ancient globular clusters host abundant and rich populations of RR Lyrae compared to the RR Lyrae-poor intermediate-age clusters, a feature easily explained by an exclusively old channel for RR Lyrae formation. Evidence that a small fraction of RR Lyrae may arise from progenitors only a few gigayears old was recently obtained from thin-disk metalrich RR Lyrae observed in the Gaia data (Zinn et al. 2019;Prudil et al. 2020;Iorio & Belokurov 2021), and from measurements of companion masses greater than 1 M e in wide-orbit RR Lyrae binaries (Kervella et al. 2019a(Kervella et al. , 2019b)). The thin-disk metal-rich RR Lyrae population was suggested by Iorio & Belokurov (2021) to be a manifestation of binary evolution pulsators (Karczmarek et al. 2017), which will register as younger stars. However, the existence of such an intermediate-age channel for LMC RR Lyrae (assuming singlestar evolution) would be in tension with the age-metallicity relation of the LMC constrained by multiple studies of field and cluster stars (Cole et al. 2005;Carrera et al. 2008;Rubele et al. 2012;Meschin et al. 2014). According to these studies, the LMC star formation history between 2 and 8 Gyr was associated with [Fe/H] between -0.4 and -1. In contrast, the LMC RR Lyrae are predominantly metal-poor, with [Fe/H]<-1 and peaking at [Fe/H]∼-1.5, as confirmed by spectroscopic studies of field RR Lyrae (Gratton et al. 2004;Borissova et al. 2006), photometric light curves of RR Lyrae (Haschke et al. 2012;Wagner-Kaiser & Sarajedini 2013;Skowron et al. 2016), and RR Lyrae-hosting globular clusters (Figure 9). Thus the LMC stellar population at ages below 8 Gyr, where we measure a significant signal in the DTD, is more metal-rich than the metallicity range measured for LMC RR Lyrae.
The DTD represents an empirical connection between the RR Lyrae sample and the SAD map. Since the RR Lyrae sample is highly complete (Section 2.1) and has strong independent evidence of originating from old stars (Section 5.1), and since our DTD recovery method would have correctly recovered a purely old RR Lyrae DTD signal from this SAD map (Section 4.1), it may be possible that the intermediate-age signal in the DTD is indicative of some systematic uncertainty in measuring older stellar ages in the LMC. The source and magnitude of this uncertainty are not obvious. The global star formation history measured by HZ09 is broadly consistent with the interaction history of the Magellanic Clouds derived from modeling of proper motion (Lin et al. 1995;Zaritsky & Harris 2004;Besla et al. 2007) and with the history of star formation and chemical enrichment derived from star clusters (Chilingarian & Asa'd 2018). As shown in Section 4.2, we verified that any incompleteness or statistical uncertainties in the MCPS photometry is unlikely to be affecting our DTD because: (1) we propagate the reported uncertainties in the SAD map into our DTDs; (2) the DTD retains signal below 8 Gyr even when measured outside the crowded LMC Bar; (3) wiping out the DTD signal of younger progenitors would require an unreasonably large unseen stellar mass with age >8 Gyr; and (4) we directly detect a signal in the DTD below 4 Gyr, where the main-sequence turnoff is detectable above the MCPS completeness limit.
On the other hand, the SAD solutions per region depend on the overall methodology adopted by the study. An example of this can be found by comparing the SAD map of the SMC from Harris & Zaritsky (2004), which was derived using roughly the same methodology used for the LMC, with the SAD map of the SMC from Rubele et al. (2018), which was derived from the deeper VISTA near-infrared survey of the Magellanic System (VMC, Kerber et al. 2009;Cioni et al. 2011). The 2-3 Gyr SAD of Harris & Zaritsky (2004) has a distinct ring pattern in the SMC, which they suspected was the result of either photometric incompleteness or systematic uncertainties in the photometric zero-point of the central SMC stars. In the VMC SAD, this ring pattern is absent in the 2-3 Gyr stellar population. Figure 11 of Rubele et al. (2018) shows that the global stellar mass formed at ages >1 Gyr in the SMC differs by more than a factor of 2 between Rubele et al. (2018) and Harris & Zaritsky (2004), and the bimodal star formation history at 2.5 and 10 Gyr found by Harris & Zaritsky (2004) is replaced with a single broad peak at 5 Gyr in Rubele et al. (2018). These systematic differences could stem from differences in initial mass function and distance assumed, the stellar isochrone models used, and/or age-metallicity binning. For example, the age-metallicity relation measured by the two studies diverges for stars older than 4 Gyr; this is most likely because Rubele et al. (2018) used metallicity bins that extend to lower abundances than Harris & Zaritsky (2004). Similar systematic differences in the LMC SADs may also exist. For example, the star formation history of the Northern Void region measured by Meschin et al. (2014), using V-and I-band photometry with the CTIO Blanco 4 m telescope, differs from that of HZ09 for ages younger than 4 Gyr. Partial estimates from upcoming SMASH data reveal a well-mixed SAD in the LMC for ages Another source of systematic uncertainty in the SADs may be the use of single-star evolution isochrone models. According to Moe et al. (2019), the close (<10 au) binary fraction for Milky Way field stars of LMC metallicity is ∼30% (compared to ∼20% for solar metallicity), so the influence of the physics of binary evolution in older populations may be non-negligible. Stanway & Eldridge (2018) have shown that models of integrated spectra and photometry of globular clusters and elliptical galaxies that include the physics of binary evolution yield age estimates that differ by a few gigayears from those of single-star models. A well-known observational manifestation of binary interaction in old stellar populations is the appearance of blue straggler stars, which are formed from the merger of ∼1 M e stars. Blue straggler stars can mimic younger stars in color-magnitude diagrams because they appear brighter and bluer than the main-sequence turnoff (Santana et al. 2016), and single-star population models that correct for the presence of blue straggler stars yield globular cluster ages older by a few gigayears (Fan & de Grijs 2012). Correcting for blue straggler stars, however, is non-trivial because their frequency is likely a function of stellar density (Santana et al. 2013(Santana et al. , 2016;;Weisz et al. 2014), and the contribution of blue straggler stars at ages of about a gigayear has been difficult to determine in composite stellar populations (Surot et al. 2019).
Leaving binary evolution aside, there are uncertainties even in the physics of single-star evolution models that can affect age estimates. For example, Tayar et al. (2017) showed that the mixing length parameter-commonly used to approximate convection theory in 1D stellar models-appears to be correlated with metallicity, and if left unaccounted for when estimating ages from the giant branch (which is the case for ages >4 Gyr in the HZ09 maps), it can lead to age uncertainties of up to a factor of 2.
Any number of these reasons could distort the SAD solutions in a subset of the cells, and therefore the final DTD, which is derived from these data.
From our work, we have shown that the DTD provides a new rigorous and quantitative test of SADs of Local Group galaxies, in addition to its original purpose as a diagnostic of stellar evolution. Although it is possible that unknown sources in systematic uncertainties in the SAD map may be driving our DTD result, we would need further tests to verify the authenticity of this issue, which we will perform in subsequent papers. Estimating the precise form of the HZ09 SAD map that would be consistent with a purely old RR Lyrae DTD is non-trivial because of the large number of SAD parameters involved in this exercise (the stellar masses per age and metallicity bin per cell), and also because RR Lyrae can only reliably constrain the oldest ages. In addition, the coarse age and metallicity binning of the HZ09 map may be less than ideal for constraining the production rate of RR Lyrae if they are coming from a narrower range of ages and metallicities within each bin (as indicated, for example, in Section 4.3. However, we can check whether our DTD result persists when calculated with SADs derived from deeper photometric data with finer age and metallicity resolution, such as the upcoming SMASH star formation histories (Ruiz-Lara et al. 2020), as well as in other Local Group galaxies such as the SMC (Rubele et al. 2018), M31 (Williams et al. 2017), and Local Group dwarfs (Weisz et al. 2014). A more constraining test of the SADs at intermediate and younger ages can also be obtained by calculating DTDs of younger variable stars with well-constrained ages, such as δ-Scutis (ages 1-3 Gyr, Petersen & Christensen-Dalsgaard 1996) and classical Cepheids (ages 70-200 Myr, Bono et al. 2005). We will pursue these in future papers.
We have calculated the first delay-time distribution (DTD) of RR Lyrae stars using the large sample of LMC RR Lyrae from the OGLE-IV survey (Soszyński et al. 2016) and the LMC's stellar age distribution (SAD) map from Harris & Zaritsky (2009). Our DTD, shown in Figure 3 and Table 1, constrains the age distribution of the full LMC RR Lyrae population, given the measured SAD of the LMC. The OGLE-IV RR Lyrae sample, which overlaps the SAD map of HZ09, contains 29,810 objects, allowing us to recover a DTD with an unprecedented balance of age resolution and detection significance. We determined the DTD signal in each age bin with an MCMC solver, and used a randomization technique to propagate uncertainties in the SAD map into the final DTD.
Our measured RR Lyrae DTD has statistically significant (>5σ) power in all age bins above 1.3 Gyr, with about 51% of the RR Lyrae associated with ages between 1.3 and 8 Gyr, and only 46% with ages above 8 Gyr (the conventional lower limit to RR Lyrae age; note that while the lower limit quoted in the literature is 10 Gyr, the SAD map has a single indivisible age bin of 8-12 Gyr, and so we refer to the lower limit as 8 Gyr in this paper). This would imply that the progenitors of RR Lyrae have zero-age main-sequence masses 2 M e at LMC metallicity, in contrast with existing constraints.
We checked the DTD for possible sources of bias. Completeness of the RR Lyrae sample is probably not an issue based on their I-band luminosity function and the predominance of RRab (fundamental) pulsators, which are least susceptible to confusion with other types of variables. We also tested our DTD algorithm on fake RR Lyrae maps drawn from a DTD that assumes all RR Lyrae are older than 8 Gyr, and found that our MCMC algorithm recovers this old DTD without any outlying detections at younger ages. The spatial distribution of RR Lyrae from a purely old DTD is also inconsistent with the spatial distribution of the OGLE-IV RR Lyrae. A possible caveat to our result is the incomplete photometry in the MCPS data and heavy crowding in the LMC central region, which limits the reliability of the SAD maps to ages younger than 4 Gyr; information about older populations is all based on HST-derived SADs in a few narrow fields. However it is not readily obvious how this is producing an intermediate-age signal since: (1) we recover the DTD signal at ages younger than 8 Gyr even after excluding the Bar region, (2) we measure a DTD signal younger than 4 Gyr, and (3) we find that in order to affect the solution, the estimates of the LMC's old stellar mass must be more than an order of magnitude greater than current measurements.
The direct interpretation of our result would be that RR Lyrae have an intermediate-age progenitor channel of 1.3-8 Gyr stars, in addition to the conventional route via ancient stars, but this possibility is in tension with existing constraints. In both the Milky Way and Magellanic Clouds, RR Lyrae are abundantly hosted in ancient globular clusters rather than intermediate-age clusters (although this conclusion is somewhat affected by smallnumber statistics due to the lower masses of intermediate-age clusters). In addition, the 1.3-8 Gyr RR Lyrae population would likely have [Fe/H]>-1 based on existing constraints on the LMC age-metallicity relation, whereas the metallicity distribution of LMC RR Lyrae measured from spectroscopy and photometric light curves is in the region with [Fe/H]<-1 with a peak at [Fe/H]∼-1.5.
The other possibility of the intermediate-age DTD result is the presence of unknown systematics in the LMC SAD map. This is not obvious because the global star formation history based on the SAD is consistent with the LMC-SMC-Milky Way interaction histories and the chemical enrichment history of the LMC derived from independent studies. However, comparison of the SMC SADs of Harris & Zaritsky (2004) and Rubele et al. (2018) shows that spatial solutions of the SAD maps can be influenced by the overall methodology adopted for their construction (e.g., assumptions about the initial mass function, the spatial size of cells, the size of age and metallicity bins, and stellar isochrone models). Age estimation of old stellar populations from color-magnitude diagrams can also be affected by binary evolution processes, such as mass transfer and mergers, as well as approximations in single-star evolution models. Any combination of these reasons could be skewing the SADs for old stellar populations, and this error could be propagating into the DTD results. We laid out further tests to investigate the physical nature of the systematics that are driving the intermediate-age signal in the DTD, such as revisiting the RR Lyrae DTD once LMC SADs from deeper photometric studies (e.g., SMASH) become available, and also measuring DTDs in the SMC and dwarf galaxies, which have less crowded star fields and deep HST-derived SADs. We will also continue to apply this technique to other types of variable stars, such as Cepheids and δ-Scutis, to probe other enigmatic phases of stellar evolution.
We are grateful to Horace Smith for reading the manuscript and providing many helpful comments and insight on RR Lyrae observations and models. We also thank the anonymous referee for many insightful comments regarding constraints on the age and chemical abundances of RR Lyrae in the Milky Way and Magellanic Clouds. We also thank Dennis Zaritsky, Knut Olsen, Thomas Matheson, Benjamin Williams, J.J. This work made use of the publicly available OGLE-IV variable star catalog.foot_3 This research has made use of NASA's Astrophysics Data System, and the VizieR catalog access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in A&AS 143, 23. This work made use of the IPython package (Pérez & Granger 2007), SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020), matplotlib, a Python library for publication quality graphics (Hunter 2007), and Astropy, a community-developed core
http://ogle.astrouw.edu.pl/
The Astrophysical Journal, 912:140 (15pp), 2021 May 10 Sarbadhicary et al.
https://people.sissa.it/~sbressan/parsec.html. We use PARSEC because it is one of the latest stellar evolution codes with both publicly available mainsequence and horizontal branch tracks. We get similar results with the MIST evolutionary tracks ofChoi et al. (2016).
http://ogledb.astrouw.edu.pl/~ogle/OCVS/