The Progenitor Systems of Classical Novae in M31

TL;DR

This paper identifies two distinct populations of classical nova progenitors in the Andromeda galaxy based on their ages and nova production rates.

Contribution

The study provides the first statistical characterization of nova progenitor systems in M31 using a delay time distribution.

Findings

A population of nova progenitors aged 2-3.2 Gyr has a nova production rate of ~3.7×10⁻⁹ events per solar mass.

Another population aged 7.9 Gyr to the age of the Universe has a higher nova production rate of ~4.8×10⁻⁹ events per solar mass.

Abstract

We present the first characterization of the statistical relationship between a
large sample of novae in M31 and their progenitor stellar populations in the
form of a delay time distribution. To this end, we leverage the spatially
resolved stellar age distribution of the M31 disk derived from deep Hubble Space
Telescope photometry by the Panchromatic Hubble Andromeda Treasury survey and a
large catalog of novae in M31. Our delay time distribution has two statistically
significant detections: one population of nova progenitors with ages between 2
and 3.2 Gyr and an unnormalized rate of (\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{upgreek}
\usepackage{mathrsfs}
\setlength{\oddsidemargin}{-69pt}
\begin{document}
${3.7}_{-3.5}^{+6.8}\pm 2.1)\times {10}^{-9}$\end{document}3.7−3.5+6.8±2.1)×10−9 events M⊙−1, and another of ages between 7.9 Gyr
and the age of the Universe with (\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{upgreek}
\usepackage{mathrsfs}
\setlength{\oddsidemargin}{-69pt}
\begin{document}
${4.8}_{-0.9}^{+1.0}\pm 0.2)\times {10}^{-9}$\end{document}4.8−0.9+1.0±0.2)×10−9 events M⊙−1 (uncertainties are statistical and
systematic, respectively). Together with the upper limits we derive at other
time bins, these detections are consistent with either a constant production
efficiency or a higher production efficiency of novae at earlier delay
times.

Figures

(a) Peak brightness in various filters vs. outburst date (used as a proxy for discovery date) for all novae in the Pietsch catalog. Novae that fall in the PHAT footprint are highlighted in pink. The average and standard deviation of the brightness in time bins of 3000 days are presented as the horizontal lines and shaded regions, respectively. (b) Time taken to decay from peak brightness to 2 mag below peak brightness vs. outburst date. (c) 300 days rolling average of the observed nova rate as a function of outburst date, compared with estimations of the global nova rate in M31 from the literature (see text for details).

The global star formation history of M31, in units of mass formed, as measured by B. F. Williams et al. (2017) and rebinned according to our temporal binning scheme. The different isochrone models yield slightly different measurements, with BaSTI emerging as an outlier in both the fifth and sixth bins. Due to their small relative size, uncertainties are omitted.

Posterior probability distributions on the DTD (in our chosen time binning scheme) for all four isochrone models. BaSTI is a notable outlier in its detection of a signal in the 600–1000 Myr bin, its lack of a detection in the 2–3.2 Gyr bin, and its disagreement on the value of the DTD in the 7.9–14.1 Gyr bin.

The PHAT SAD map of M31, calculated using MIST and binned according to our temporal scheme, with the Pietsch nova catalog overplotted. A darkly shaded spatial cell indicates higher star formation rate, with the shading normalized separately for each time bin.

Comparisons between the SAD map derived from the four isochrone models for selected time bins. The gray scale is normalized independently for each plot.

(a) Violin plot of our posterior distributions for the DTD rate in each time bin. (b) Violin plot of our completeness-corrected DTDs plotted against theoretical DTDs from A. J. Kemp et al. (2021) and an observed SN Ia DTD from D. Maoz et al. (2012). The DTD in blue is derived using the effective survey length; the DTD in purple is derived using the time-limited subset of the nova sample.

Full text

Introduction

Novae are among the most common astrophysical transients arising from binary stellar evolution. However, no large-scale study of their progenitors has yet been conducted. A key concern in the study of binary systems is their effect on the formation and evolution of compact objects: in binaries with different initial masses, one star will always leave the main sequence first, giving rise to interactions between a living star and a stellar remnant. For example, when a star begins to transfer matter from its outer layers onto its degenerate companion—whether by overfilling its Roche lobe or through strong stellar winds (P. Podsiadlowski & S. Mohamed 2007)—it creates the conditions for a thermonuclear explosion on the surface of the white dwarf. These explosions eject the envelope of accreted mass but are nonterminal for the remnant and the companion star. They often produce an optical transient, commonly known as a classical nova (S. Starrfield et al. 1972; J. S. Gallagher & S. Starrfield 1978; D. Prialnik & A. Kovetz 1995; D. M. Townsley & L. Bildsten 2004—see L. Chomiuk et al. 2021, for a recent review). The details of the mass transfer process, the build-up to the thermonuclear runaway, and the consequences they both have for the subsequent evolution of the components of the binary system are complex. Inquests into this process, whether by stellar evolution codes (P. A. Denissenkov et al. 2012, 2014; B. Paxton et al. 2015) or binary population synthesis (BPS) simulations (H.-L. Chen et al. 2016; A. J. Kemp et al. 2021), are only possible by making considerable simplifying assumptions that significantly cloud our picture of the landscape of nova progenitors.

Given their relative abundance compared to other astrophysical transients, novae provide a highly accessible probe into binary stellar evolution, enabling the use of a rich statistical toolset on larger samples to constrain the nature of their progenitor systems. A precise measurement of the evolutionary timescales and formation efficiencies of nova progenitors would provide observational constraints to help answer a number of open questions in binary stellar evolution. Among these questions are the influence of the initial conditions (since it is now well established that the fraction of close binary systems in the main sequence is a strong function of stellar properties such as mass and metallicity; see M. Moe & R. Di Stefano 2017; C. Badenes et al. 2018; M. Moe et al. 2019; C. N. Mazzola et al. 2020), the impact of the orbital and stellar parameters on the stability of mass transfer (K. Pavlovskii et al. 2017; K. D. Temmink et al. 2023), the role played by stellar winds in the mass transfer process (S. Mohamed & P. Podsiadlowski 2007; N. A. Webb 2023), and the many uncertainties involved in the onset, progression, and aftermath of common-envelope episodes (N. Ivanova 2011). Many of these uncertainties are encoded explicitly or partially in BPS codes, which have been used to study novae by H.-L. Chen et al. (2016, 2019) and A. J. Kemp et al. (2021, 2022). A common thread in these theoretical studies has been a dearth of observational constraints derived from large, statistically significant samples of novae.

To address these issues, we present here the first measurement of the delay time distribution (DTD) of novae in M31. The DTD is the occurrence rate of a class of objects as a function of time following a single brief burst of star formation. By characterizing the spatial correlation (or lack thereof) between the ages of stars and the objects of interest, we can recover the formation efficiency of that class of objects as a function of lookback time in field stellar populations—another way of defining the DTD (D. Maoz & C. Badenes 2010; C. Badenes et al. 2015; S. K. Sarbadhicary et al. 2021). An observationally derived DTD can be used to test theoretical expectations, including the predictions from BPS models and their underlying assumptions.

The Andromeda galaxy (M31), the closest large galaxy to the Milky Way, is the ideal environment for a large-scale study of nova progenitors. At a distance of 752 ± 27 kpc (A. G. Riess et al. 2012), the stellar populations in M31 can be resolved by the Hubble Space Telescope (HST) down to magnitude ∼27 in regions of low stellar density, which has allowed the Panchromatic Hubble Andromeda Treasury (PHAT) team (J. J. Dalcanton et al. 2012; B. F. Williams et al. 2014, 2017) to use precise multiband photometry to produce spatially resolved maps of the stellar age distribution (SAD) in an area that encompasses roughly one third of the M31 disk. We have combined the PHAT data with the extensive historical catalog of novae in M31 from W. Pietsch et al. (2007), which contains over a thousand entries, to derive the DTD for novae in M31. This paper is structured as follows: Section 2 reviews the nova catalog and PHAT SAD map in detail, Section 3 explains the theory behind the DTD and the process of recovering it, Section 4 presents our DTD in the context of the stellar isochrones used to generate the SAD map, and Section 5 explores our results in more depth and compares them to previous literature.

Data

Our work builds on decades of observations of M31; namely, a historical nova catalog covering the entire galaxy, originally created for comparison with an X-ray nova catalog (W. Pietsch et al. 2007), but regularly updated on a publicly available website maintained by W. Pietsch.7 ^7^ https://www.mpe.mpg.de/m31novae/opt/m31/ We combine these data with a spatially resolved ancient star formation history of M31 (B. F. Williams et al. 2017) produced from the PHAT survey (J. J. Dalcanton et al. 2012; B. F. Williams et al. 2014). This map provides the crucial link between the locations of novae and the age distribution of the stars surrounding them, enabling the recovery of a DTD.

Nova Catalog

2.1.

The nova catalog goes as far back as the seminal survey of Andromeda led by Edwin Hubble (E. P. Hubble 1929) and other surveys conducted using photographic plates (N. U. Mayall 1931; W. Baade & H. H. Swope 1963; C. Payne Gaposchkin 1963). More recent additions come mainly from the Zwicky Transient Facility (E. C. Bellm et al. 2019), which scans the entire optical northern sky every 2 days. This high-cadence survey has a depth of 20.6 in the R band. The catalog is also populated by an array of observations from independent astronomers searching for transients and reporting their discoveries on databases such as the Central Bureau for Astronomical Telegrams and the Transient Name Server. By their nature, these disparate additions to the nova inventory are recorded in a variety of filters and limiting magnitudes. As of 2022 November, the catalog contains 1219 individual entries, 263 of which fall in the footprint covered by the PHAT survey (see Section 2.2 below).

Because the nova catalog is heterogeneous and spans more than a century of observations, we must examine its contents critically. It is likely that some novae might have been missed in early studies or that faint and/or fast novae are underrepresented in surveys conducted before the advent of high-efficiency CCD detectors and high-cadence surveys, leading to biases and completeness issues. To evaluate these issues, we show the peak brightness and decay times of the novae in the Pietsch catalog as a function of discovery date in Figure 1, as well as a rolling average of the observed nova rate. The peak brightness of M31 novae ranges from magnitude 13.9 to 20.7, with an average and standard deviation of 17.3 ± 1.0. The decay times—only reported for a subset of the sample—span a much wider range between 2 and 410 days; taking the average and standard deviation in log space yields 1.28 ± 0.42 log(days).

(a) Peak brightness in various filters vs. outburst date (used as a proxy for discovery date) for all novae in the Pietsch catalog. Novae that fall in the PHAT footprint are highlighted in pink. The average and standard deviation of the brightness in time bins of 3000 days are presented as the horizontal lines and shaded regions, respectively. (b) Time taken to decay from peak brightness to 2 mag below peak brightness vs. outburst date. (c) 300 days rolling average of the observed nova rate as a function of outburst date, compared with estimations of the global nova rate in M31 from the literature (see text for details).

To look for evidence of evolution in either the peak brightness or decay timescales as a function of discovery date, we divide the novae into temporal bins of 3000 days for the former and 6000 days for the latter (see Figure 1). Our analysis shows no statistical evidence for such evolution.

The measured nova rate has grown significantly in recent years, reflecting the increased efficiency of nova surveys. In Figure 1(c), we compare the rolling average of the nova rate to three estimates from the literature: those of M. Capaccioli et al. (1989) and A. W. Shafter & B. K. Irby (2001), which focused on the historical record of observed novae in M31, and that of M. J. Darnley et al. (2006), who corrected for completeness using artificial nova tests and models of M31’s surface brightness and internal extinction in the context of the POINT-AGAPE survey (S. Calchi Novati et al. 2005). The rolling average of the observed nova rate in the most recent epochs included in the W. Pietsch et al. nova catalog (∼​​​​​​30 novae yr^−1^) is consistent with the values derived by M. Capaccioli et al. (1989) and A. W. Shafter & B. K. Irby (2001), but a factor 2 lower than the value of M. J. Darnley et al. (2006). This latter discrepancy is to be expected, as the rate reported in that study accounts for novae that would be unobservable and it cannot be directly compared to observed rates. We conclude that while the W. Pietsch et al. nova catalog is certainly incomplete, it does not appear to have significant biases against any subclass of novae included in the sample, such as faint or fast-evolving novae.

Another property of the catalog entries relevant to our DTD analysis is the spatial accuracy of the nova positions. In theory, uncertainties on nova sky coordinates could propagate to uncertainties on nova counts in each spatial cell (see Section 2.2). Unsurprisingly, novae discovered at earlier times tend to have larger uncertainties. The earliest members of the nova catalog were discovered on photographic plates, with their positions (and accompanying uncertainties) being reported as a Cartesian distance from some chosen center as opposed to sky coordinates (A. W. Shafter et al. 2015). However, the median uncertainty of the entire catalog is half an arcsecond—corresponding to 1.8 pc (projected) at the distance of M31, far smaller than the 300 × 1400 pc (deprojected) spatial cells of the SAD map—and no nova has a coordinate uncertainty greater than 13″, or 47 pc (projected). We conclude that errors in the spatial location of M31 novae that could lead to “cell hopping” and propagate into uncertainties on the DTD are likely rare.

To remove contamination by recurrent novae (i.e., counting separate outbursts of the same system as different systems), we rely on A. W. Shafter et al. (2015), who screened the W. Pietsch et al. catalog looking for spatially coincident nova candidates in the original digital images or photographic plates. By removing all instances of recurrence reported by their analysis except for the most recent nova candidate, we ensure that each progenitor system gives rise to exactly one nova and every nova included in our DTD analysis is unique. Since this study was conducted, an additional 255 novae have been added to the W. Pietsch et al. catalog; to remove recurrent novae from these entries, we looked for those that showed spatial overlap within 1σ of another nova. In total, our recurrence analysis removed 31 nova outbursts from the catalog, 10 of which fall within the PHAT footprint, leaving a total of 253 unique novae to be compared with the spatially resolved SAD. We also remedied an apparent typo in the outburst date of nova 2013-02a, recorded as Julian Date 3456332.98 (January 7th, 4751 CE), changing the leading digit from 3 to 2.

In order to report nova rates as well as abundances, we must account for the clear incompleteness of the catalog at early times. We attempt this correction in two ways: by deriving an effective survey length and by only analyzing a sample of recent novae, where the catalog is consistent with measured nova rates.

To calculate an effective survey length, we assume that the catalog is effectively complete in the final time bin of 3000 days (see Figure 1(a)), which corresponds to a rate of 30 novae yr^−1^—in good agreement with other calculations of observed M31 nova rates (M. Capaccioli et al. 1989; A. W. Shafter & B. K. Irby 2001). We correct for incompleteness at earlier times with the following procedure: for each time bin of 3000 days, we calculate the ratio of the nova count in that period and the count in the final period, multiply that ratio by the time span of the bin, and increment the total effective survey length by that diminished value rather than the full 3000 days. This calculation, which depends on the uncertain choice of a “complete” nova rate against which we compare our historical sample, yields an effective survey length of 38 yr. We stress that this number is merely a best estimate—as are the nova rates derived from it—and heterogeneous historical catalogs such as this cannot be precisely corrected for completeness.

The most recent 6000 days of the W. Pietsch et al. catalog have an observed rolling nova rate consistent with previous studies (see Figure 1). By limiting our sample to these recent novae, we can provide a rough cross-check to our effective correction for survey length completeness. However, this approach shrinks our nova sample considerably and limits our ability to recover the DTD. The results of both completeness analyses are presented in Section 5.

PHAT Stellar Age Distribution Map

2.2.

The PHAT survey collected HST photometry for 117 million individual stars in the northern half of M31, measuring the stars simultaneously in six bands with coverage from the UV to the IR. This rich data set, hosted on MAST at doi:10.17909/T91S30 (J. Dalcanton & B. Williams 2012), has resulted in dozens of publications that focus on different aspects of the stellar content of M31. B. F. Williams et al. (2017) recovered the star formation history by dividing the PHAT data into 826 83″ × 83″ (deprojected to 0.3 × 1.4 kpc at the distance of M31) spatial cells, producing color–magnitude diagrams (CMDs) for each cell—assuming a Kroupa initial mass function (P. Kroupa 2002)—that were then fit with the well-tested MATCH software package (A. E. Dolphin 2002). MATCH also contains tools to model both systematic (A. E. Dolphin 2012) and random (A. E. Dolphin 2013) uncertainties, which the authors implement when deriving their SAD map.

Certain features in a CMD are clear indicators of the presence of stars of a certain age; by comparing an observed CMD to a model composed of simple stellar populations of various ages, an SAD can be robustly recovered. B. F. Williams et al. (2017) used MATCH 2.6 to calculate the best-fit SADs to CMDs drawn from several bands, only including stars brighter than the 50% completeness limit and making use of the software’s ability to incorporate complex extinction models. They draw from a high-resolution dust map of M31 (J. J. Dalcanton et al. 2015) derived from a two-component model of reddening-induced spread among red giant branch stars. This model provides a galactic map of median AV, spread of the log-normal AV distribution, and fraction of stars being affected by dust; B. F. Williams et al. (2017) then feed these parameters to MATCH. They also account for foreground reddening by fitting an additional parameter AV FG independently for each spatial cell.

Each cell was fit using four different sets of isochrones: Padova (P. Marigo et al. 2008; L. Girardi et al. 2010), BaSTI (A. Pietrinferni et al. 2004, 2013; S. Cassisi et al. 2006), PARSEC (A. Bressan et al. 2012), and MIST (J. Choi et al. 2016). This resulted in four maps of the SAD (or star formation history) in the PHAT footprint, one for each set of isochrones, each with 16 temporal bins spanning from the present day to the age of the Universe (see B. F. Williams et al. 2017 for details). The metallicities used to generate the isochrones were determined not by fitting a single Z value to each spatial cell, but by modeling the enrichment history of the entire galaxy. Upon finding evidence for chemical enrichment at early times (M. Mollá et al. 1997), the authors make three radial divisions to M31 and fit an exponentially decreasing (in time) enrichment rate to each region. Their best-fit enrichment decay rates are independently consistent with regions of high star formation and the observed metallicity gradient in M31.

We assigned each of the 253 unique novae in the W. Pietsch et al. catalog that fall within the PHAT footprint to one of the 826 spatial cells in the B. F. Williams et al. (2017) maps. The median PHAT cell has no novae, but 158 cells—about 19%—have at least one, with the two most populated cells containing nine distinct novae each. To find a good compromise between retaining enough temporal resolution in the DTD and maximizing the likelihood of detecting power in at least some of the combined bins, given the size of our sample, we combined the 16 time bins of the native PHAT SAD maps into seven larger bins. We retained the youngest (0–300 Myr) and oldest (7.9–14.1 Gyr) temporal bins. The remaining native bins in the SAD maps, which vary in duration from 100 Myr to >1 Gyr, were merged so as to be roughly equal in logarithmic time. This rebinning scheme improved the average ratio between Mi,j, the stellar mass formed in spatial cell i and time bin j, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\sigma }{{M}{i,j}}\end{document} (the uncertainty on that quantity), while greatly reducing the number of spatial cells with Mi,j = 0. The SAD obtained with this temporal binning scheme, integrated over the entire PHAT footprint for each of the four isochrone sets, is shown in Figure 2.

The global star formation history of M31, in units of mass formed, as measured by B. F. Williams et al. (2017) and rebinned according to our temporal binning scheme. The different isochrone models yield slightly different measurements, with BaSTI emerging as an outlier in both the fifth and sixth bins. Due to their small relative size, uncertainties are omitted.

The four isochrone data sets used by B. F. Williams et al. (2017), and the corresponding SAD maps, offer a unique opportunity to evaluate the impact that systematic uncertainties associated with stellar evolution models have on the star formation histories recovered from resolved stellar populations. This topic is treacherous and uncertain (see C. Gallart et al. 2005; C. Conroy et al. 2009, for discussions), but it remains a fundamental limitation in attempts to draw statistical inferences about the association between resolved stellar populations and the products of stellar evolution, as we will show in Section 5.

Methods

In recovering the DTD, this work follows the procedure laid out by D. Maoz & C. Badenes (2010) and C. Badenes et al. (2015). Convolving the DTD (Ψ(t), with units of events yr^−1^ M⊙^−1^) with the SAD (M(t), with units of M⊙) yields an event rate over time R(t) (units of events yr^−1^):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}R(t)={\int }_{0}^{t}M(\tau ){\mathrm{\Psi }}(t-\tau )d\tau .\end{eqnarray*}\end{document}

In practice, we only have access to the current event rate R(t0) and an SAD split into discrete time bins; however, the PHAT survey and M31 nova catalog give us this information for 826 spatial cells. We can then recast the integral in Equation (1) discretely, as a dot product:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{R}_{i}({t}_{0})=\displaystyle \sum _{j=0}^{{t}_{0}}{M}_{i,j}{{\mathrm{\Psi }}}_{j}\,\end{eqnarray*}\end{document}

where Ri(t0) is the modern-day nova rate in the ith spatial cell, Mi,j is the total stellar mass formed in the ith cell and jth temporal bin, and Ψ_j_ is the DTD, or event rate per year per unit mass, in the jth temporal bin (which is inherent to novae and therefore common across all spatial cells). The DTD recovery can now be treated as an inverse problem, where the values of the DTD in each temporal bin, Ψ_j, are free parameters that are varied to achieve the set of rates Ri(t0) that best fit the observed nova counts ni_ in all bins, simultaneously.

Typically, Ri(t0) would be calculated by dividing the number of novae in each spatial cell, ni, by the effective survey length (D. Maoz & C. Badenes 2010). As noted in Section 2.1, we estimate an effective survey length of 38 yr for the W. Pietsch et al. catalog. However, due to the uncertainty in this estimate, we adopt ni as a proxy for Ri(t0) and report our DTD in units of events M⊙^−1^, nova count per unit of stellar mass, rather than nova rate per unit stellar mass (except where otherwise noted).

To explore this seven-dimensional parameter space and obtain Bayesian posteriors on each Ψ_j_, we employ the dynamic nested sampling routine dynesty (J. S. Speagle 2020; S. Koposov et al. 2024). In short, nested sampling (J. Skilling 2004, 2006) is an efficient method of recovering posteriors P(θ∣Data, Model) on a set of parameters θ by simultaneously estimating the Bayesian evidence P(Data∣Model) and the posterior. Dynamic nested sampling (DNS; E. Higson et al. 2019) allocates the posterior samples adaptively throughout the sampling process, exploring high-likelihood regions of the parameter space more thoroughly. DNS has many attributes that make it an attractive alternative to methods such as the Markov Chain Monte Carlo process, chief among them being its speed and the trivial independence of its samples.

In our application of DNS, we employ a log-likelihood function drawn from the modified chi-squared statistic derived by K. J. Mighell (2020):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{\chi }_{\gamma }^{2}=\displaystyle \sum _{i=0}^{N}\frac{{({n}_{i}+{\mathrm{\min }}({n}_{i},1)-{m}_{i})}^{2}}{{n}_{i}+1}\,\end{eqnarray*}\end{document}

where ni is the observed number of novae in a spatial cell and mi is the expected number of novae according to our model. The latter number is calculated by convolving the model DTD with the observed SAD in the ith spatial cell. We implement the Bayesian routine with the single bounding method of P. Mukherjee et al. (2006) and the slice sampling method of R. M. Neal (2003) and W. J. Handley et al. (2015a, 2015b).

The modified chi-squared statistic, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\chi }{\gamma }^{2}\end{document} , is designed to converge to the true mean value of Poisson-distributed data (even in the low-N regime), unlike the standard χ^2^. We employ this modified statistic rather than drawing a log-likelihood function directly from the Poisson distribution in order to directly incorporate variance in the expected number of novae mi_ deriving from the statistical variance in the SAD used to calculate it. Instead of propagating this uncertainty through bootstrap or hybrid Monte Carlo methods (A. E. Dolphin 2013), we exploit the fact that the denominator of the summed term in Equation (3) corresponds to the variance of ni and add the variance on mi, which is calculated through uncertainty propagation ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\mathrm{Var}}({m}{i})={\sum }{j=0}^{{t}{0}}{\sigma }{{M}{i,j}}^{2}{{\mathrm{\Psi }}}{j}^{2}\end{document} ).

We then minimize \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\chi }_{\gamma }^{2}\end{document} by maximizing the log-likelihood of the corresponding Gaussian distribution,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\displaystyle \begin{array}{rcl}{\mathrm{ln}}(L) & = & -\frac{1}{2}\displaystyle {\sum }_{i=0}^{N}\frac{{\left({n}_{i}+{\mathrm{\min }}({n}_{i},1)-\displaystyle {\sum }_{j=0}^{{t}_{0}}{\dot{M}}_{i,j}{{\mathrm{\Psi }}}_{j}\right)}^{2}}{{n}_{i}+1+\displaystyle {\sum }_{j=0}^{{t}_{0}}{\sigma }_{{M}_{i,j}}^{2}{{\mathrm{\Psi }}}_{j}^{2}}\\ & & +\,{\mathrm{ln}}\left(2\pi \left({n}_{i}+1+\displaystyle {\sum }_{j=0}^{{t}_{0}}{\sigma }_{{M}_{i,j}}^{2}{{\mathrm{\Psi }}}_{j}^{2}\right)\right),\end{array}\end{eqnarray*}\end{document}

applying DNS to this equation, and finding the parameters Ψ_j_ that best fit our data. We note that this “χ^2^-to-Gaussian” approach is also taken by B. D. Johnson et al. (2021) in their spectral energy distribution-fitting code prospector.

Any Bayesian routine requires the establishment of a prior distribution on the parameters being measured. We employ the same prior on all Ψ_j: flat between 0 and 1.5 × 10^−6^ events M⊙_^−1^ and zero elsewhere. This upper limit is the rate necessary to produce the observed 1σ upper limit on M31 nova rate (M. J. Darnley et al. 2006) in the case that the DTD is only nonzero in the single time bin with the least star formation.

For our final results, we calculated the edges of the 5% highest probability density (HPD) region for each parameter and took the midpoint as the maximally likely rate. Similarly, our 1σ uncertainty ranges on each parameter correspond to the 68.27% HPD region, following C. Badenes et al. (2015). If the lower limit of this latter region coincides with the lower edge of the parameter space (technically, if the region falls within the smallest 100 parameter values explored by the dynesty fitting routine, with a typical DNS run involving upwards of 10^5^ such values), we report the rate to be statistically consistent with zero—a nondetection—and instead record a 2σ upper limit.

Results

In Figure 3, we present the posteriors on event rate per unit stellar mass in each time bin, computed separately for the four available isochrone models. The DTD is a mathematical representation of the spatial correlation between events and stellar populations, as illustrated in Figure 4. In this figure, produced using the MIST isochrone set, there is a clear correlation between nova positions and stellar mass of ages 2–3.2 Gyr and 7.9–14.1, which indicates the possibility of a physical link between those stars and that progenitor population of novae. We expect to see—and indeed do see—significant detections in the DTD corresponding to these temporal bins. We also note that our DTD is dissimilar to a DTD generated from random sky positions, showing that these detections are specific to novae.

Posterior probability distributions on the DTD (in our chosen time binning scheme) for all four isochrone models. BaSTI is a notable outlier in its detection of a signal in the 600–1000 Myr bin, its lack of a detection in the 2–3.2 Gyr bin, and its disagreement on the value of the DTD in the 7.9–14.1 Gyr bin.

The PHAT SAD map of M31, calculated using MIST and binned according to our temporal scheme, with the Pietsch nova catalog overplotted. A darkly shaded spatial cell indicates higher star formation rate, with the shading normalized separately for each time bin.

The most prominent feature of these posteriors is the robust recovery of a population of nova progenitors in the oldest time bin, with ages between 7.9 Gyr and the age of the Universe, though BaSTI yields a detection that is ∼40% lower than the average of the other three. This divergence is unsurprising, as B. F. Williams et al. (2014) found BaSTI to be most discrepant from the other models at very old ages. According to MIST stellar evolution models, these lifetimes correspond to zero-age main sequence (ZAMS) masses of 1.06–0.90 M⊙. With the exception of BaSTI, the other three isochrone sets (MIST, Padova, and PARSEC) also yield a detection in the 2–3.2 Gyr time bin. BaSTI and PARSEC also detect a signal in both the 0.6–1 Gyr and 1–2 Gyr time bins, but MIST and Padova do not. Disagreements at these young ages are not unexpected; the models become less accurate as the main-sequence turnoff point falls below the photometric depth of the PHAT survey. For all other time bins, the maximally likely event rate was either less than the standard deviation of the posterior (which we treat as a nondetection) or essentially zero.

There is broad agreement between the isochrone sets, and the results are consistent with our expectation of higher formation efficiencies at earlier delay times from BPS models (discussed further in Section 5). We note that the exact values, not just the general features, of the DTDs from each isochrone set are broadly consistent when dynesty is run multiple times on the same SAD map. For our final DTD measurement, we take statistical and systematic uncertainties into account by averaging only the results of the PARSEC, Padova, and MIST posteriors and combining the standard deviations of the three posteriors in quadrature, then presenting these uncertainties alongside the standard deviations between the three models (see Table 1). Having excluded BaSTI from the final DTD, we can see that the statistical uncertainties are larger than the systematic differences between isochrone sets, although the range of variation is significant: the ratio of statistical to systematic uncertainty is ∼2:1 in the 2.0–3.2 Gyr bin, but ∼9:1 in the 7.9–14.1 Gyr bin. The final row of Table 1 shows that, as anticipated, the expected number of novae in all spatial cells derived from our recovered DTD, ∑imi, falls within 1σ of the observed number of unique novae in the PHAT footprint ∑ini, which is 253.

Table 1: Maximum Likelihood and 1σ Uncertainty Ranges of the Posteriors for the DTD from Each Model (Units of Events M⊙−1)

The most notable disagreements between the DTDs derived from the four isochrone sets, in the time bins where they do differ, can be traced back to disagreements on the location of star formation in specific time bins. An illustrative example of this is the 600–1000 Myr bin shown in Figure 3. The signal grows more statistically significant in the order MIST, Padova, PARSEC, BaSTI—corresponding to the sequence of increasing stellar mass formed in that same time bin (see Figure 2).

Figure 5 further clarifies this disagreement. The concentration of star formation in the 10 and 20 kpc rings, which is present in all four models in the first time bin (0–300 Myr), yields a consistent nondetection in the corresponding bin of Figure 3, given that the spatial distribution of novae shows no enhancement at these locations. However, the varying amount of star formation in the 20 kpc ring recovered by each isochrone set in the 600–1000 Myr range leads to the noted discrepancies in the recovered DTD. In the 1–2 Gyr time bin, the MIST and Padova models detect substantial star formation in the outer regions of the disk, whereas PARSEC and BaSTI only detect it close to the bulge—colocated with the vast majority of novae. For this reason, the latter two models report a statistically significant nova rate in the corresponding time bin of the DTD, and the former two do not. These discrepancies between the DTDs derived using the different isochrone sets stress the model-dependent nature of DTD analyses and the importance of taking into account systemic biases and uncertainties.

Comparisons between the SAD map derived from the four isochrone models for selected time bins. The gray scale is normalized independently for each plot.

In Figure 6, we present a violin plot of our final DTD (in events M⊙^−1^) and two completeness-corrected DTDs: one obtained by dividing our results by the effective survey length of 38 yr (obtained in Section 2.1), and one obtained by limiting our sample to recent novae. These two DTDs are statistically consistent with each other except in the final time bin, where the 1σ bounds come close to overlapping, but do not. We compare these completeness-corrected DTDs with a set of theoretical DTDs from A. J. Kemp et al. (2021), as well as an observationally determined DTD for Type Ia supernovae (SNe Ia) from D. Maoz et al. (2012), which we discuss in the following section.

(a) Violin plot of our posterior distributions for the DTD rate in each time bin. (b) Violin plot of our completeness-corrected DTDs plotted against theoretical DTDs from A. J. Kemp et al. (2021) and an observed SN Ia DTD from D. Maoz et al. (2012). The DTD in blue is derived using the effective survey length; the DTD in purple is derived using the time-limited subset of the nova sample.

Discussion and Conclusions

We have presented the first statistical inference of the formation efficiency of novae across stellar populations of different ages (i.e., a DTD for novae) in M31. Our results are derived using SAD maps in M31 from deep HST photometry produced by the PHAT collaboration (B. F. Williams et al. 2017), and the catalog of M31 novae compiled by W. Pietsch, which spans more than 100 years of observations. Taking statistical and systematic uncertainties into account, our DTD has robust detections in two time bins: stars between 2 and 3.2 Gyr of age and between 7.9 Gyr and the age of the Universe. While our ability to recover signal in certain time bins may be due to the particular nature of M31’s star formation history and nova catalog, the DTD itself is intrinsic to the event of interest and, critically, does not depend on the data used to recover it.

The fact that M31’s global SAD is dominated by older stars (see Figure 2) naturally makes the recovery of a DTD signal at young ages more difficult, particularly with a sample of modest size. Additionally, the excess of stellar mass in the 2–3.2 Gyr age range facilitates a detection at these delay times; one possible explanation for this signal is a burst in nova progenitor formation driven by the 2–3 Gyr old galactic merger (F. Hammer et al. 2018). In this context, the question of which stellar populations generate most of the novae in M31 is distinct from the question of which stellar populations are most efficient at generating novae in general.

Taken together, our detections and upper limits are consistent with a rate of nova formation that decays with delay time, although, strictly speaking, we cannot rule out a uniform formation rate. The maximally likely rate of formation efficiency we recover is, perhaps surprisingly, higher in the final time bin than in the 2–3 Gyr bin, possibly implying a departure from monotonic decay of nova formation. In fact, there is no theoretical expectation that the DTD of novae should be exactly monotonic—the theoretical DTD calculated by A. J. Kemp et al. (2021) exhibits several local maxima. However, given the wide range of rates falling within our 1σ uncertainty intervals, we caution against overinterpretation of this result.

A decaying formation rate is borne out by simple arguments in stellar evolution. Technically, a binary system can produce novae as soon as the more massive partner evolves into a white dwarf (WD). As delay time increases, the ZAMS masses of both the donor and the degenerate accretor decrease on average, leading to a decrease in formation efficiency of novae (H. Ritter et al. 1991; U. Kolb 1995; H.-L. Chen et al. 2016) and an increase in recurrence time (O. Yaron et al. 2005), as less mass is available to transfer to the accretor. Although the PHAT SAD map lacks sufficient resolution at short delay times to probe the “turn-on” phase—the prediction that the formation efficiency of novae should be zero before the formation of the first WDs at delays of ∼40 Myr—it does probe, and successfully recover, the broad expectation of a decrease in formation efficiency at longer delays.

Our rough completeness corrections allow us to estimate nova rates (per unit stellar mass) as well as nova counts; these rates are consistent with the shape and height of the A. J. Kemp et al. (2021) DTD (see Figure 6), with the effective survey length correction yielding a slightly better match than the time-limited sample. To be clear, no adjustments have been made to our data to bring them into alignment with the theoretical rates. This agreement is simply a product of the DTD recovery process outlined in Section 3 and the completeness corrections presented in Section 2.1. In the same figure, we also make a comparison against the Kemp et al. DTD of helium novae—rare events caused either by the deposition of helium, rather than hydrogen, onto a WD or by the fusion of hydrogen into helium on the surface of the dwarf before a nova detonation. As expected, He nova rates are so low as to be inconsistent with our completeness-corrected DTD; this large discrepancy implies that a minute fraction of novae in M31 are He novae.

Previous BPS studies of novae, such as by H.-L. Chen et al. (2016) and A. J. Kemp et al. (2021), provide theoretical support for more efficient nova production at earlier times. The claim made by Chen et al., that less massive binaries will generate fewer novae, is borne out by their nova rate over time for two model SFHs. A. J. Kemp et al. (2021) go into much greater depth with regard to the progenitor population, providing detailed descriptions of common evolutionary pathways and their implications for nova rates. They present a (downward-sloping) DTD of their own and split it up by the evolutionary phase of the donor star. Following from Figure 12 of A. J. Kemp et al. (2021), the progenitor accretor population of the 2–3.2 Gyr time bin should be dominated by O/Ne WDs. The progenitor donor population is a more eclectic mix, with stars from the low-mass main sequence and first giant branch best represented. The accretor population in the 7.9–14.1 Gyr bin should be an approximately even mix of O/Ne and C/O WDs, with first-ascent giant branch stars dominating the donor population.

Studies of novae, both in our Galaxy and in Andromeda, have been plagued by dust-driven spatial incompleteness, the extent of which has been a subject of ongoing debate. M31’s bulge-dominated nova population stands in contrast to that of other galaxies (R. Ciardullo et al. 1989) and the previously discussed expectation that younger stellar populations should be more efficient progenitors of novae. Were our catalog to be biased against novae in the disk, any correlation between young stellar populations and novae would be confounded, artificially driving down the DTD at early delay times. K. Hatano et al. (1997) used a simple model of dust in M31 to argue that the observed bulge-to-disk nova ratio is a consequence of dust extinction in the disk, and the true fraction of novae residing in the bulge is of the order of 25%. However, their conclusion is deeply model-dependent; even their limited explorations of changes to this model are consistent with a bulge-dominated nova population. Notably, a map of extinction in M31 created from the PHAT survey (J. J. Dalcanton et al. 2015) does not resemble the Hatano et al. model.

A. W. Shafter & B. K. Irby (2001) use the spatial distribution of planetary nebulae in M31, which should be slightly more sensitive to extinction than novae, as a tracer of dust. They find this population to be less centrally concentrated than novae, concluding that the bulge dominance of the nova population in M31 is genuine and not an artifact of the dusty disk. There is no doubt that the W. Pietsch et al. catalog misses some novae due to dust. However, the result from Shafter & Irby and the fact that our results are consistent with higher formation efficiency among younger stellar populations make the existence of a large missing population of disk novae that would introduce severe spatial incompleteness unlikely.

Recent literature suggests a population of “faint, fast” novae (M. M. Shara et al. 2017) that would evade surveys with cadences much longer than an hour. This prediction, also supported by the theoretical models of O. Yaron et al. (2005), has not been borne out in newer surveys. Between 2018 and 2019, the Zwicky Transient Facility surveyed fields in the Galactic plane with a typical cadence of 40 s to search for short-period astrophysical variables (T. Kupfer et al. 2021); no such novae have been reported. Shara et al. predict that, at the distance of M31, these objects would have magnitude 17–18 and decay times t2 of 5 hr or less. Our nova catalog certainly probes these magnitudes, but appears not to be sensitive to such short delay times (see Figure 1). Therefore, our DTD does not contain information about these transients, which may not even comprise a significant fraction of novae.

Ultimately, one can only ever produce a DTD for observed events. The kinds of light-curve and extinction completeness corrections undertaken in studies of absolute nova rates can estimate the total number of novae left unobserved, but naturally cannot estimate the location of each unobserved nova—the data that would be required to recover a DTD.

We conclude with a brief discussion of the comparison between our nova DTD and the DTD for Type Ia SNe recovered by D. Maoz et al. (2012), which is interesting because some subclasses of novae have been proposed as potential SN Ia progenitors (D. Maoz et al. 2014). Our completeness-corrected nova DTD is about four orders of magnitude higher than the SN Ia DTD, and noticeably shallower at long delay times. This implies a strong (≲0.1%) upper limit on the fraction of nova-producing systems that go on to explode as SNe Ia in star-forming galaxies like M31.

Such a discrepancy presents a problem for the single-degenerate theory of Type Ia progenitors, which posits that a significant fraction of SNe Ia arise from gradual mass transfer between a living star and a remnant—the same kind of system that produces novae. If such a small fraction of these systems end their lives as type Ia supernovae (right panel of Figure 6), there may not be enough single-degenerate systems in star-forming galaxies to explain the observed rate of SNe Ia (D. Maoz et al. 2014).

Our work reinforces the previously demonstrated power of DTDs for constraining the progenitor populations of the products of binary stellar evolution. The continuation of high-cadence surveys of nearby, well-studied galaxies such as M31 would provide better constraints on the nova DTD at all lookback times and reduce our reliance on historical nova catalogs of dubious completeness. Resolved stellar populations of those galaxies provide opportunities for the recovery of DTDs for other astrophysical transients. The upcoming extension of the PHAT SAD map to the southern region of M31 (Z. Chen et al. 2025) promises a near-trivial doubling of the nova sample size; once published, these data will represent an exciting opportunity to improve the resolution of the nova DTD.

SAD maps similar to that of B. F. Williams et al. (2017) exist for other galaxies in the Local Group (e.g., M. Lazzarini et al. 2022 for M33, J. Harris & D. Zaritsky 2009; A. Mazzi et al. 2021 for the LMC, J. Harris & D. Zaritsky 2004; S. Rubele et al. 2018 for the SMC), with matching, high-quality catalogs of astrophysical transients and other products of stellar evolution. This research enables further application of the methods described here, which could be leveraged to constrain BPS models, deepening our understanding of key stages of binary stellar evolution and improving the predictions for the rates of rare events such as SNe Ia and black hole mergers.

The data and code underlying this work will be shared upon a reasonable request made to the lead author.