wdwarfdate: A Python Package to Derive Bayesian Ages of White Dwarfs

White dwarfs have been successfully used as cosmochronometers in the literature, however their reach has been limited in comparison to their potential. We present wdwarfdate, a publicly available Python package to derive the Bayesian age of a white dwarf, based on its effective temperature (Teff) and surface gravity (logg). We make this software easy to use with the goal of transforming the usage of white dwarfs as cosmochronometers into an accessible tool. The code estimates the mass and cooling age of the white dwarf, as well as the mass and main-sequence age of the progenitor star, allowing for a determination of the total age of the object. We test the reliability of the method by estimating the parameters of white dwarfs from previous studies, and find agreement with the literature within measurement errors. By analyzing the limitation of the code we find a typical uncertainty of 10% on the total age when both input parameters have uncertainties of 1%, and an uncertainty of 25% on the total age when Teff has an uncertainty of 10% and logg of 1%. Furthermore, wdwarfdate assumes single star evolution, and can be applied to calculate the total age of a white dwarf with parameters in the range 1,500<Teff<90,000 K and 7.9<logg<9.3. Finally, the code assumes a uniform mixture of C/O in the core and single star evolution, which is reliable in the range of white dwarf masses 0.45-1.1 Msun (7.73<logg<8.8).


INTRODUCTION
White dwarfs are the perfect candidates to use as cosmochronometers. These degenerate objects are formed with high temperatures from the death of stars with masses < 8 M , and cool over time. The cooling process is understood to the point that by deriving the effective temperature (T eff [K]) and surface gravity (log g [cm s −2 ], units for log g will be implicit from now on to simplify notation) from observational data we can estimate how long it has been cooling with a relatively high accuracy (e.g., Fontaine et al. 2001;Bédard et al. 2020). These cooling models, combined with stellar evo-lutionary models (e.g., Choi et al. 2016;Dotter 2016) and the initial-to-final mass relation (IFMR, e.g., Cummings et al. 2018;Marigo et al. 2020), provide a reliable method to estimate the total age of a white dwarf. In addition, with the release of Gaia DR2 and EDR3 (Gaia Collaboration et al. 2016Brown et al. 2021), a large number of new white dwarfs have been discovered -about 10 times more than previously known-that have precise proper motions and parallaxes measured by Gaia (e.g., Gentile Fusillo et al. 2019McCleery et al. 2020).
The Gaia mission also allowed for the identification of binaries with unprecedented precision (e.g., El-Badry et al. 2021). Given that it can be assumed that the components of a binary system were born at approximately the same time, identifying white dwarfs with wide co-moving companions allows us to calculate the age of the system under the assumption that each object has likely evolved as a single star, without phases of mass transfer or common envelope evolution (Bodenheimer 2011). Stellar age is one of the most difficult fundamental properties to measure directly; although some methods such as gyrochronology (e.g., Skumanich 1972;Barnes 2003;Van Saders et al. 2016), asteroseismology (e.g., Chaplin et al. 2014;Aguirre et al. 2017) and isochrone fitting (e.g., Baraffe et al. 2015;Agüeros et al. 2018;Berger et al. 2020) can yield relatively precise ages for single stars, they are limited to specific ranges of masses, or stages of mainsequence evolution, in which they can be applied (e.g., Chabrier & Baraffe 1997; Baraffe et al. 2015;Rodríguez et al. 2016). For example, none of these methods can be applied to low-mass stars (< 0.6 M ). Binary stars provide a powerful tool to estimate stellar ages "independently of mass": the age of a star of any mass can be estimated if it is co-moving with another object which can be age-dated.
White dwarf cosmochronology is the subject of extensive literature. Some examples include Kalirai (2012), who calculated the ages of white dwarfs in the Milky Way halo to study its formation; Gagné et al. (2018), who estimated the age of the ultra-massive white dwarf GD 50 to refine the age of the AB Dor moving group which GD 50 belongs to; Anguiano et al. (2017), who estimated white dwarf parameters, including their ages, to study the formation of the Galactic disk; Catalán (2015) who used white dwarf-white dwarf binaries to calibrate the low-mass end of the IFMR; and Kilic et al. (2019), who used white dwarf ages to study the Milky way inner halo. In addition to these, other studies (e.g., Garcés et al. 2011 andFouesneau et al. 2019) have estimated the age of main-sequence stars by calculating the age of co-moving white dwarf companions. Furthermore, previous studies have presented software libraries that can be used to estimate white dwarf ages: WD models 1 (S. Cheng, priv. comm.) is a Python library to estimate white dwarf parameters from Gaia photometry, although it does not provide measurement uncertainties for these parameters; and BASE-9 2 (von Hippel et al. 2006(von Hippel et al. , 2014) is a C ++ library that estimates star cluster and stellar parameters from photometry, and in particular can be used to estimate white dwarf and progenitor masses and ages. BASE-9 uses a Bayesian technique, combined with a Markov chain Monte Carlo and numerical integration techniques to estimate the posterior probability distributions of the different parameters, including age. As a consequence, BASE-9 results complicated to implement, in particular when only the parameters of the white dwarf are needed. This shows the incredible potential of estimating total ages of white dwarfs, and the lack of tools to make it accessible.
In this paper, we present wdwarfdate, a new Python library to estimate the total age of a white dwarf with the associated uncertainty, from T eff (K), log g, and atmospheric composition. The goal of this software is to make the usage of white dwarf as cosmochronometers accessible. Therefore, this code is designed to be easily used, with warnings for the corresponding limitations of the theoretical models applied in the process. The estimation of the white dwarf age is done using Bayesian statistics under the assumptions of a white dwarf core with uniform mixture of C/O and single star evolution. wdwarfdate also calculates the mass and cooling age of the white dwarf and the mass and main-sequence age of its progenitor star, along with their respective uncertainties (and their full probability densities) which include scatter in the IFMR and propagate the measurement errors of the white dwarf parameters used as an input. This paper is divided into the following sections: in Section 2, we describe the general method used to estimate the total age of a white dwarf (Section 2.1), the models that were included for the age determination (i.e., cooling tracks, IFMRs and stellar evolutionary models; Section 2.2), the Bayesian framework on which wdwarfdate relies (Section 2.3), and the fast-test mode, an alternative method to estimate the parameters of the white dwarf included in the code (Section 2.4). In Section 3, we apply wdwarfdate to calculate the total ages of a reference sample of age-calibrated white dwarfs to validate our method. We also estimate the ages of 18 white dwarfs in wide co-moving pair candidates with M dwarfs (Section 3.3) to provide a sample of age-calibrated M dwarfs for future studies. In Section 4, we detail the constraints and assumptions that are required in using wdwarfdate. Our results are summarized in Section 5.

METHOD
In this section, we describe the sequential steps required to calculate the mass and cooling age of the white dwarf, as well as the mass and the main-sequence age of the progenitor star, and finally the total age of the object. wdwarfdate will derive ages using Bayesian statistics -Bayesian age-by default, but the code also includes a fast-test mode, an alternative method to estimate the parameters of a white dwarf, and we explain both in this section. We also describe the models in which wdwarfdate relies for the parameter estimation. For an example on how to use wdwarfdate see Appendix A.

Outline of the algorithm
The calculation of the total age using wdwarfdate relies on a precise measurement of T eff and log g. The estimation of these two parameters is done from observational data, and it is described in Section 2.2. From T eff and log g, the process to estimate the total age of a white dwarf is divided into four stages, summarized in Figure 1: 1. From the T eff and log g input values combined with cooling evolutionary tracks for the white dwarf we obtain its mass (final mass, m f ) and cooling age (t cool ).

2.
From the m f we estimate the mass of the progenitor star (initial mass, m i ) using an IFMR, which are, in most cases, empirically calibrated.
3. With the m i combined with stellar evolution tracks we estimate its evolutionary lifetime, meaning the time since the star formed until it became a white dwarf, which we will call main sequence age or lifetime in this work (t ms ).
4. Finally, adding t cool and t ms , we obtain the total age of the object (t tot ).
These steps depend on models which assume single star evolution. If a white dwarf interacted with another object at any point, then the evolutionary models will not reproduce its parameters correctly (e.g., Marsh 1995;Brown et al. 2011). In Section 4 we discuss different scenarios where the single star evolution assumption might not be valid. Steps required to estimate a white dwarf total age from an effective temperature and surface gravity (T eff and log g). In this diagram, m f and mi are the final mass (mass of the white dwarf) and initial mass (mass of the progenitor star), and t cool , tms and ttot are the cooling age of the white dwarf, the main sequence age of the progenitor star and the total age of the object, respectively. σx indicates the uncertainty of each parameter x.

Models incorporated in wdwarfdate
The accuracy and precision of the estimated fundamental white dwarf parameters using wdwarfdate depend strongly on how the two input observables were determined, the obtained values, and uncertainties. In addition, to initialize wdwarfdate the user needs to choose between cooling models for DA white dwarf and non-DA, an IFMR and a stellar evolution model for the progenitor star (metallicity and rotation). In this section, we describe the methods generally used to estimate T eff and log g in the literature and how they affect predictions by wdwarfdate. We also describe the models included in wdwarfdate for each of the choices described above, as summarized in Table 1.
The determination of temperature and surface gravity generally rely on either a photometric or a spectroscopic approach. Photometric methods rely on accurate trigonometric parallaxes and high precision photometry (e.g. Koester et al. 1979) to reconstruct the spectral energy distribution of a white dwarf, and requires an accurate knowledge of the extinction caused by interstellar dust. Spectroscopic methods, however, rely on atmosphere models to reproduce the detailed profiles of hydrogen or helium lines, and depend on their theoretical line profiles (e.g. Bergeron et al. 1992;Beauchamp et al. 1999). These line profiles in most cases approximate the white dwarf atmosphere as a 1-dimensional object, which can lead to biases in some extreme cases of white dwarf properties, and in most cases ignore the impact of magnetic fields which can be significant for the most massive white dwarfs. Recent studies using data from Gaia DR2 and the spectroscopic method show that the results from these two techniques generally agree within the uncertainties , except for T eff 6500 K where the parallax-based photometric estimations of log g were found to be more accurate than spectroscopic estimations (Napiwotzki et al. 2020;Kawka & Vennes 2012).
wdwarfdate relies on the theoretical evolutionary sequences of the Montreal White Dwarf Group . These models assume that the white dwarf has a core composed of a uniform carbon and oxygen mixture (X C = X O = 0.5), surrounded by a He layer (with a fractional mass q He = M He /M * = 10 −2 ) surrounded by an outermost H layer. The evolutionary tracks assume that this outermost H layer is "thick" (q H ≡ M H /M * = 10 −4 ) for hydrogen-atmosphere white dwarfs (DA white dwarf), and "thin" (q H = 10 −10 ) for helium-atmosphere white dwarfs (non-DA white dwarf; Bédard et al. 2020). The publicly available Montreal cooling tracks were calculated for masses in the range 0.2−1.3 M in steps of 0.05 M . These tracks are shown in Figure 2 in black, together with the limits within which wdwarfdate was designed to function in terms of T eff and log g, shown in a blue dashed line. These limits are based on the following limitations: (1) the lowest temperature covered by the cooling tracks (1500 K); (2) the full range of white dwarf masses covered by the cooling tracks (0.2 − 1.3 M ); and (3) the 0.3 Myr isochrone, given that cooling ages below this value are unreliable. In summary, the approximate limits within which wdwarfdate will rely on the cooling tracks are: 1, 500 T eff 100, 000 K and 7.0 log g 9.3. Section 4 details additional limitations of wdwarfdate.
To estimate the progenitor mass from the white dwarf mass, or vice versa, the user can choose among the following IFMRs: a spectroscopic determination of the white dwarfs T eff and log g, which was then used to estimate the cooling age and mass of the white dwarf using the Fontaine et al. (2001) cooling tracks. These models are valid for T eff 30, 000 K, and were updated by the more recent Bédard et al. (2020) cooling tracks. By subtracting the cooling age from the age of the host open cluster, they obtained the main-sequence age of the progenitor star, which they then used to estimate the progenitor mass using the PAdova and TRieste Stellar Evolution Code (PARSEC, Bressan et al. 2012) and the models computed using the Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018: the MESA Isochrones & Stellar Tracks (MIST, Choi et al. 2016;Dotter 2016). These two families of stellar evolutionary models resulted in two distinct IFMRs, both of which are included in wdwarfdate. The difference between the initial masses calculated for the two IFMR are generally within 5%. By default, wdwarfdate uses the MIST-based IFMR as recommended by Cummings et al. (2018), due to the nonrotating MIST isochrones estimating more precise ages for young stars (∼ 100 Myr), as tested with the Pleiades cluster. However, the PARSEC-based IFMR should be considered especially for high-mass progenitor stars (> 5 M ), given that the MIST models tend to underestimate masses in this range.
The semi-empirical IFMR from Marigo et al. (2020), was also obtained from white dwarfs from known clusters. In particular, they included seven new white dwarf members of the old open clusters NGC 752 (1.55 Gyr) and Ruprecht 147 (2.5 Gyr), which allowed them to study the low-mass end (around 0.5 M ) of the IFMR in detail. In addition, Marigo et al. (2020) used a new analysis technique which combines photometric and spectroscopic data to estimate precise T eff , log g and m f , and to confirm cluster membership for the white dwarfs. As shown in Figure 3, they found a kink at low-masses, which they associated with the formation of carbon stars in the Galaxy. The IFMR from Salaris et al. (2009) also used white dwarfs from known clusters to calibrate the relation, but put special emphasis in the models used, and performed a detailed analysis of the sources of uncertainties and scatter in the relation. Finally, the IFMR from Williams et al. (2009) studied white dwarfs from the intermediate-age open cluster M35 (NGC 2168), which provided information to constrain the higher-mass end of the IFMR.
The IFMRs included in wdwarfdate depend on cluster objects that are younger and therefore more massive than average field white dwarfs, which could cause a bias in the relation. In addition, as some of the IFMRs were done prior to the Gaia mission, and many of these objects were assumed to be cluster members based on their kinematics and colours, without parallaxes, some of these objects might not be true members, biasing the relation. The IFMR will also place an extra boundary in the limits of T eff and log g for which wdwarfdate can estimate a total age. For example, Cummings et al. (2018) MIST-based IFMR was calibrated for white dwarf masses between 0.5 − 1.3 M , which will place a tight constraint specially in the lowest value possible of log g. We will discuss these limitations further in Section 4.
The IFMRs from Cummings et al. (2018) and Marigo et al. (2020) are calibrated for m i > 0.83 M , and we included them in the code extending this limit to 0.45 M , assuming a constant m f between 0.45 < m i < 0.83 (gray line in Figure 3). This extension allows wdwarfdate to estimate parameters for white dwarfs with masses m f 0.6 M , which is the peak of the white dwarf mass distribution (e.g., Genest-Beaulieu & Bergeron 2019).
To estimate the main sequence age from the initial mass, or progenitor mass, we adopted in wdwarfdate the MIST isochrones. These tracks cover a mass range between 0.1 − 300 M with a step of 0.05 M , and an age range from 0.1 Myr to 13.8 Gyr. To select the optimal stellar evolution track, the user has to choose the metallicity and rotation of the progenitor star among the following options: [Fe/H] = {−4, −1, 0, 0.5} (dex), [α/Fe] = 0, and v/vcrit = {0.0, 0.4}. We will discuss the effects of using different tracks on the total age in Section 4.    In this section, we describe the Bayesian method implemented in wdwarfdate. Our code derives Bayesian ages of white dwarfs, with a posterior probability given by

Bayesian framework
where m i is the initial mass, or the mass of the progenitor star, and t cool is the cooling age of the white dwarf. ∆ m is an extra parameter which takes into account the scatter in the IFMR of white dwarfs, and it is described below. T eff and log g are the inputs for the effective temperature and surface gravity, respectively. To simplify notation we did not include uncertainties in Equation 1, although these are also required in the code. Using Bayes theorem, the posterior probability in Equation 1 can be written as where p(log 10 t cool ), p(m i ) and p(∆ m ) are the priors on the cooling age, the initial mass and ∆ m respectively, which are described below, and is the likelihood. T eff and log g in Equation 3 are the modeled effective temperature and surface gravity of the white dwarf. These modeled quantities are calculated from the three parameters that are being sampled (m i , t cool and ∆ m ) as shown by the probabilistic graphical model in Figure 4, which represents the posterior probability of the problem. In summary: 1) From the initial mass (m i ) we estimate a main sequence age (t ms ); 2) by adding t ms to the cooling age (t cool ) we obtain the total age (t tot ) of the white dwarf; 3) from m i and the delta parameter (∆ m ) we estimate a final mass (m f , described below); 4) from m f and t cool we estimate the modeled effective temperature and surface gravity (T eff and log g) which are compared to the input parameters (T eff and log g) to estimate the posterior probability. The models used in each step are described in Section 2.2. The ∆ m parameter in Equations 1, 2 and 3 was included to model the scatter in the IFMR. This scatter has different causes including contamination from noncluster members, measurement errors, metallicity variations, inconsistencies with the models used and age determination of the clusters (e.g., Weidemann 2000; Williams et al. 2009;Casewell et al. 2009). The parameter ∆ m is drawn from a normal distribution centered at zero with standard deviation σ m , so that We estimated the parameter σ m from the Cummings et al. (2018) data, and adopted σ m = 0.03 M . The inclusion of the scatter in the posterior is indicated in Figure 4 asm f = m f + ∆ m , where m f is the true mass of the white dwarf andm f is the scattered mass, because it was calculated using the IFMR. The priors indicated in the posterior probability of Equation 2 are given by the Star Formation History (SFH) for t cool and the Initial Mass Function (IMF) for m i . We assume a constant SFH (e.g., Madau & Dickinson 2014) and applied it as a flat prior on t tot , where we required t tot < 15 Gyr to avoid a sharp cut at the age of the Universe in the age probability distributions. As t tot = t cool + t ms , the Jacobian of the coordinate transformation is 1, and this prior could be used without extra modifications. In addition, we assumed the IMF to be m −2.3 i (Kroupa 2001;Chabrier 2003), which describes the range of stellar masses expected to form a white dwarf. In short, T eff and log g are calculated from the parameters being sampled (mi, t cool and ∆m), and these values are compared to the input parameters. In the figure,m f is the scattered final mass (mass of the white dwarf) calculated from the initial mass using an IFMR, and ∆m is an extra parameter to take into account the scatter in this relation. Therefore, m f is the true final mass, calculated by adding ∆m tom f .
wdwarfdate obtains the probability distributions by explicitly integrating over the likelihood and priors to perform the normalization. We set the code to evaluate the posterior in a grid made of discrete arrays for each parameter: m i , log 10 t cool and ∆ m . The code first performs a wide grid evaluation with the parameters ranging between 0.3 Myr < t cool < 15 Gyr, 0.1 M < m i < 10 M and −0.1 M < ∆ m < 0.1 M . Then it obtains the probability distribution for each parameter by setting the minimum and maximum evaluation limits according to where the probability is higher. These limits can also be adjusted by the user. We use these results combined with the models described in Section 2.2 to derive probability distributions for the rest of the parameters (total age, final mass and main sequence age).

Fast-test method
One of the limitations of wdwarfdate is given by the IFMR: the total age of a white dwarf can be estimated only in the cases where its mass is between the allowed limits of the IFMR. For the cases where a total age cannot be estimated but the T eff and log g are within the limits of the cooling tracks, we included in wdwarfdate a fast-test method to estimate the final mass and cooling age only. For this method, wdwarfdate generates two normal distributions for the effective temperature and surface gravity, using the input values (T eff and log g) as means, and the uncertainties (σ T eff and σ log g ) as the standard deviations of the distributions such as The code does a Monte Carlo uncertainty propagation by running the two distributions through the process described in Section 2.1 and in Figure 1, and obtains a distribution for each of the parameters of interest: final mass, initial mass, cooling age, main sequence age and total age. wdwarfdate will automatically run this method when the final mass of the white dwarf is outside of the allowed ranges by the IFMR. The main difference between the fast-test and the Bayesian methods are the priors: the fast-test method does not include the IMF and SFH as priors on the initial mass and the total age, respectively, but it does constrain the total age to be smaller than 15 Gyr. To test that the final mass and cooling age estimated with the fast-test and the Bayesian methods are comparable, we run a grid of white dwarfs with T eff = 1, 500 − 100, 000 K and log g = 7 − 9.3 and uncertainties of 10% and 1% respectively, using the two methods. For both methods we used the DA white dwarf cooling sequences, the Cummings et al. (2018) MIST-based IFMR, the MIST tracks with [Fe/H] = 0 and v/vcrit = 0.0, and the uncertainties were calculated using the 16 th and 84 th percentile of the distributions. The results show that the final mass estimations agree within the uncertainties between both methods (right panel Figure 5), and the cooling age estimations agree in most cases (left panel Figure 5). We color-coded in purple the points for which the cooling age differs by more than 1σ, and we found that these points agree in the final mass within the uncertainty, although the cooling age values seem to be more discrepant as the cooling age decreases. This discrepancy does not depend on the cooling age because there are some white dwarfs with small cooling ages which have the same result in the fast-test and Bayesian method. The purple points on Figure 5 have some of the highest temperatures and surface gravities, as shown by the (T eff − log g) inset diagram in the left panel of Figure 5. Moreover, the purple points follow the area where the cooling tracks are closer to each other in Figure 2. This increases the uncertainty in the cooling age, making the likelihood probability distribution function (PDF) of these objects wider. Therefore, the prior, in particular the IMF, has a bigger effect on the estimated values, making the cooling ages estimated with the Bayesian method larger, given that it favors smaller initial masses.
The final masses between 0.65 − 0.8 M in the right panel of Figure 5 are slightly lower with the Bayesian method than with the fast-test method. This effect is also caused by the prior, in a similar way as described above. By inspecting the cooling tracks in Figure 2, we noticed that the cooling tracks in the mass range 0.65 − 0.8 M are closer to each other, making the likelihood less constrained. The differences for the final mass and cooling age described above show the advantage of applying the Bayesian method which adds extra information using the priors (IMF and SFH), unlike the fasttest method which provides the highest likelihood value.
In summary, when the uncertainties in T eff and log g are small (∼ 10% and ∼ 1%, respectively) the fast-test and the Bayesian method are comparable in the estimation of final masses and cooling ages, except for white dwarfs with both T eff > 40, 000 K and log g > 8.4, which coincide with dense track areas. However, these exceptions represent only 2% of our sample, and some have extremely high temperatures, therefore these cases are not common. Fast-test method Figure 5. Comparison of the fast-test and Bayesian methods included in wdwarfdate for the cooling age and final mass (left and right panels, respectively). We included the oneto-one relation in a black solid line. We run a grid of white dwarfs with T eff = 1, 500 − 100, 000 K and log g = 7 − 9.3 and uncertainties of 10% and 1%, respectively. This comparison shows that the estimation of final mass and cooling age using both methods is similar. The purple points in both panels represent calculations of cooling age that differ more than 1σ between the two methods. The small panel of T eff and log g shows that these points are situated in a dense area of isochrones (Figure 2), which makes the likelihood less constrained. As a consequence, the prior has more influence and makes the Bayesian cooling ages larger.
In this section we compare the results from wdwarfdate with those from previous studies. The goal of this section is to test the validity of wdwarfdate results. We also calculate the ages of white dwarfs in candidate binaries with M dwarfs which make a new set of age calibrators. Unless indicated otherwise, the estimations were made using the Bayesian method and the Cummings et al. (2018) MIST-based IFMR, with the DA white dwarf cooling sequences, and [Fe/H] = 0 and v/vcrit = 0.0 stellar evolution models. The uncertainties were calculated using the 16 th and 84 th percentile of the distributions.

Comparison to white dwarfs from clusters from
Cummings et al. (2018) Cummings et al. (2018) used a sample of 79 white dwarfs from 13 different clusters of known ages with T eff and log g measurements from the literature to calibrate the IFMR. In addition, the authors re-measured T eff and log g for 8 of the white dwarfs. Cummings et al. (2018) is one of the most comprehensive studies to-date of white dwarfs from clusters. They used the cooling tracks from the Montreal White Dwarf Group to estimate final masses and cooling ages from T eff and log g, and the MIST isochrones to estimate initial masses from the cooling and total ages of the objects (See Section 2.2 for a detailed description of their work). Given that they estimated all the parameters which can be obtained with our code using a similar procedure, Cummings et al. (2018) is an ideal study to compare the performance of wdwarfdate.
We run wdwarfdate on all the stars in Cummings et al. (2018) using the T eff and log g and uncertainties published in their work. We compared final mass, cooling age, initial mass and main sequence age in Figure 6. Most of the values agree within 1σ with their results. The uncertainties calculated with wdwarfdate are comparable to the estimated by Cummings et al. (2018) for the cooling age and final mass. For the main sequence age and initial mass, the uncertainties are smaller in the Cummings et al. (2018) results because they were propagated from the age of the cluster which has a small uncertainty, and not the IFMR as in wdwarfdate.
We color-coded the white dwarfs for which any of the parameters in Figure 6 did not agree within 1σ to study them in detail. The black dot is the white dwarf GD 50, and although the parameters agree within the uncertainties it is color-coded because it is discussed below. The white dwarfs color-coded in green and purple are outliers in the IFMR (Figure 3), which explains why the initial mass and main sequence age do not agree. We refer to Cummings et al. (2018) for the analysis of the outliers. The final mass and cooling age for the green objects agree within uncertainties, which is expected because we are using the same T eff and log g to estimate them. In the case of the three white dwarfs color-coded in purple, wdwarfdate estimated a higher cooling age than expected according to the parameters in Cummings et al. (2018). However, these objects have T eff > 50, 000 K, and the cooling tracks used in Cummings et al. (2018) (Fontaine et al. 2001) are not suitable for temperatures > 30, 000 K. We confirmed our calculations of the cooling ages with the Montreal White Dwarfs Data Base (MWDD, Dufour et al. 2017) 3 , which contains the updated models by Bédard et al. (2020), which we implemented in wdwarfdate.
To compare the total age estimated with wdwarfdate with the cluster age, we excluded the outliers (points color-coded in green and purple) in Figure 6. The total estimated age agrees within 1σ with the cluster age for 65% of the white dwarfs, and within 2σ for 94%. We also calculated a median relative difference of 12% between the median age for each white dwarf that was not an outlier and the age according to the cluster membership, Both the green and purple points are outliers in the IFMR. In addition, the three purple white dwarfs with higher cooling age have high effective temperatures, which are better modeled with the new cooling tracks implemented in our code . The black dot is the massive white dwarf GD 50.
which shows wdwarfdate is doing a good job estimating these total ages. The white dwarf color-coded in black is GD 50, which is a young and ultramassive white dwarf of 1.28 ± 0.08 M . For this object we obtained a main sequence age of 46.21 +4.76 −3.74 Myr, a cooling age of 82.62 +11.19 −14.08 Myr, a total age of 128.71 +10.54 −11.92 Myr, a final mass of 1.27 ± 0.02 M , and an initial mass of 7.5 ± 0.3 M . Our progenitor mass calculation differs slightly from the value in Cummings et al. (2018), but agrees with the results of Gagné et al. (2018), that calculated a mass of 7.8 ± 0.6 M . They found GD 50 is likely member of the AB Doradus moving group, and performed a detailed estimation of its parameters using the Montreal white dwarf evolutionary models and the MIST models, accounting for all possible C/O/Ne core compositions, and finding a total age of 117±26 Myr. Our result is between this value and the one estimated by Cummings et al. (2018) that assumed GD 50 belongs to the the Pleiades (based on Dobbie et al. 2006) and assigned it an age of  Figure 7. We compare the total age posterior PDF obtained with wdwarfdate for the white dwarfs in Cummings et al. (2018), with the age of the cluster they belong to. Also, we divided the white dwarfs according to their membership and added in a black vertical line the age of the group. We found good agreement between the age of the cluster and the age derived with wdwarfdate.
135 ± 35 Myr, which is slightly larger than the current measured age of this cluster (112 ± 5 Myr; Dahm 2015). We also used wdwarfdate on GD 50 using the PARSECbased IFMR from Cummings et al. (2018), as was suggested in Section 2.2 for progenitor masses > 5 M , and we obtained a total age of 114.61 +12.72 −11.68 Myr, significantly closer to the results of Gagné et al. (2018).

Comparison to white dwarfs from M67 from Canton et al. (2021)
To test the accuracy of wdwarfdate at lower masses, we estimated the parameters of the white dwarfs from Canton et al. (2021). The authors measured T eff and log g from high resolution spectra for 22 white dwarfs from the M67 cluster which is 3.5 Gyr and solar metallicity. They used the cooling tracks from Fontaine et al. (2001) to estimate final masses from these parameters, and the PARSEC isochrones to estimate initial masses from the age of the cluster. They estimated initial masses only for a smaller vetted sample, and discarded the rest because, for example, they were potential Hecore white dwarfs or potential blue straggler remnants. Using their measured T eff and log g, and uncertainties with wdwarfdate, we found good agreement for final mass and cooling age between our results and Canton et al. (2021) -as shown in Figure 8-for all but two white dwarfs. However, our results for these two white dwarfs agree with the MWDD, therefore the differences are likely due to Canton et al. (2021) using an older version of the Montreal cooling tracks (Fontaine et al. 2001). In addition, we compared our estimation of the total age and the initial mass for the vetted sample from Canton et al. (2021), and found that our results for both parameters agree within 1σ, as shown by the bottom panels in Figure 8. The uncertainties on the total age of the white dwarfs are large (±5 Gyr) which is expected given that these are low-mass white dwarfs and the largest contribution to the total age comes from the main sequence age, which tends to have higher uncertainties. Furthermore, the uncertainties on the initial mass calculated by Canton et al. (2021) are smaller because they used the known total age to estimate these masses, which is well constrained. To further check our code, we run wdwarfdate using the IFMR from Marigo et al. (2020) which differs from Cummings et al. (2018) for low-masses, and found the same results for the median value of the parameters, although the shape of the posterior was different in some cases, due to the change in the shape of the IFMR.

Wide, co-moving white dwarf companions
We used wdwarfdate to estimate ages of 30 white dwarfs identified to be candidate wide co-movers with M dwarfs by Kiman et al. (2021). These M dwarfs constitute a new set of age-calibrators and we refer the reader to Kiman et al. (2021) for information about the main-sequence co-moving star. Kiman et al. (2021) estimated the white dwarf ages using wdwarfdate, and the T eff and log g calculated by Gentile Fusillo et al. (2019) with photometry from Gaia DR2. For our calculations we used the updated version of that work (Gentile Fusillo et al. 2021), that took advantage of the precise photometry and parallaxes from Gaia EDR3 and estimated T eff and log g for 359, 073 high-confidence white dwarf candidates. However, we did not find improvement between Gaia DR2 and EDR3 in the relative uncertainty of T eff and log g for the 30 white dwarfs, and we found a relative difference for the values of 8.4% and 2.2% respectively.
All of the 30 objects have a white dwarf probability > 95% assigned by Gentile Fusillo et al. (2021). However, one of the white dwarfs does not have an estimation for T eff and log g. We discarded that source along with the white dwarfs with masses < 0.5 M , or > 1.1 M , according to their position in the color-magnitude diagram ( Figure 9) and white dwarfs with G BP −G RP > 0.9. In short, the reasons for these cuts are that high-and low-mass white dwarfs are thought to be the result of binary evolution, and therefore do not follow the single  Canton et al. (2021). These white dwarfs belong to the M67 cluster (3.5 Gyr). For the final mass, cooling age and initial mass we show the oneto-one line in black, and for the total age we show the age distributions with the age of the cluster in a black line. For the total age and initial mass we only show the vetted white dwarfs from Canton et al. (2021). We found good agreement between our results and their study. The uncertainties in the initial mass calculated by Canton et al. (2021) are smaller than the point in some cases because they used a well constrained total age to estimate these masses. star evolution assumption made by the cooling tracks included in wdwarfdate . In addition, the color-cut avoids contamination from nonwhite dwarfs such as unresolved binaries. For a more detailed discussion of the limitations of wdwarfdate see Section 4.
We used wdwarfdate to estimate the ages of the remaining 18 white dwarfs. Five of these objects were identified as DA white dwarfs in the MWDD, two as DB, and we assumed DA for the rest, and used the corresponding T eff and log g from Gentile Fusillo et al. (2021). Given that most white dwarfs are DAs (∼ 75%), this is a good assumption (Fontaine et al. 2001;Gentile Fusillo et al. 2015). To run our code we used the cooling sequence with thick H outermost layer for DA white dwarfs and with thin layer for non-DA, the stellar evolution models for . We show as red empty circles the white dwarfs which were discarded because they are outside the ranges which are likely to have followed single star evolution: 0.5 M < m f < 1.1 M , and also white dwarfs with GBP − GRP < 0.9 to avoid contamination from non-single white dwarfs. The only object within these limits which is also discarded did not have an estimation of T eff and log g in Gentile Fusillo et al. (2021). We also show the models of constant mass 0.5 M , 0.7 M , 0.9 M and 1 M for DA white dwarfs (Bergeron et al. 1995;Fontaine et al. 2001;Holberg & Bergeron 2006;Kowalski & Saumon 2006;Bergeron et al. 2011;Tremblay et al. 2011;Blouin et al. 2018).
calculated these white dwarf ages using the Marigo et al.
(2020) IFMR and found similar results, which are not displayed here. The Gaia source id of the white dwarfs, the input T eff and log g, and the parameters estimated with wdwarfdate are in Table 2. The reported values in this table are the median, 16 th , and 84 th percentile as uncertainties. We were able to estimate the total age for 13 of the 18 white dwarfs; the rest have masses outside of the allowed range by the IFMR. For these 13 white dwarfs we used both the Bayesian and fast-test methods included in wdwarfdate, and both results are shown in the table, while for the remaining 5 we only used the fast-test method to estimate the final mass and cooling age. The white dwarfs from Table 2 exemplify the impact of the T eff and log g uncertainties in the estimation of the parameters, and show the advantage of using the Bayesian method. For example, Gaia EDR3 1323951092358466304 has a 15% and 7% uncertainty on T eff and log g, respectively. Although the parameters estimated for this white dwarf are within the uncertainties between the fast-test and the Bayesian method, there is a significant difference between, for example, the final masses (0.64 +0.16 −0.06 M using the Bayesian method and 0.97 +0.23 −0.26 M using the fast-test method). In this case where the likelihood is not well constrained due to the uncertainties, the prior in the Bayesian method (in particular the IMF) is adding extra information preferring lower-mass white dwarfs, while the fast-test method selects the highest likelihood value. This was discussed in detailed in Section 2.4, and will be discussed further in Section 4.

LIMITATIONS AND ASSUMPTIONS
White dwarf cosmochronology is a powerful tool to estimate stellar ages. However, it is based on a number of models and assumptions that need to be taken into account. In this section, we discuss the constraints and model choices when using wdwarfdate.
We run a grid of 100 white dwarfs with 1, 500 K < T eff < 100, 000 K and 7 < log g < 9.3 with different relative uncertainties to test the resulting uncertainty of the Bayesian total ages derived with wdwarfdate. In this study we used the Cummings et al. (2018) MISTbased IFMR, with the DA white dwarf cooling sequences and [Fe/H] = 0 and v/vcrit = 0 stellar evolution models. We found that for uncertainties of 1% in both T eff and log g, the median relative uncertainty in total age is ∼ 11%, as shown in Figure 10. If the uncertainty of T eff is increased to 10%, the total age uncertainty increases to ∼ 25%. On the other hand, using an uncertainty of 1% for T eff and 5% for log g, the total age uncertainty increases to ∼ 52%, showing the importance of precise input measurements, in particular of log g. For a study of the total age relative uncertainty as a function of T eff and log g see Appendix B. We repeated this analysis for non-DA cooling models and found no significant difference in the results.
In addition to the limits on T eff and log g set by the cooling tracks discussed in Section 2.2, the IFMR imposes extra limits to the input values for which wdwarfdate will be able to estimate a total age, main sequence age and initial mass. To quantify this limit, we run wdwarfdate on a grid of T eff and log g and recorded which parameters were estimated. In this exercise, we adopted the MIST-based Cummings et al. (2018) IFMR, the cooling sequences for DA white dwarfs, and the stellar evolution models for [Fe/H] = 0 and v/vcrit = 0.0. The values of T eff and log g for which wdwarfdate can estimate a total age are color-coded by final mass in Figure 11. In summary, wdwarfdate will derive Bayesian total ages for 1, 500 K T eff 90, 000 K and 7.9 log g 9.3, and only final mass and cooling age using the fast-test method for 7.0 log g 7.9 (green-blue five point stars in Figure 11). The masses for C/O core white dwarfs ranges from 0.45−1.1 M . For lower masses, the core is composed of He and for higher masses of O and/or Ne (Fontaine et al. 2001;Garcia-berro et al. 1997). As we implemented cooling tracks that assume a C/O core in wdwarfdate, parameter estimation using our code outside this range should be taken with caution. We discuss further the lower and upper mass limits below.
Within a Hubble time, it is not possible to have a white dwarf with a mass smaller than 0.45 M that was formed by single stellar evolution. Low-mass white dwarfs (< 0.45 M ) cannot ignite He, therefore they have a He core, and are thought to be the result of binary evolution (e.g., Marsh 1995;Brown et al. 2011). There are studies which propose that these objects could also be formed from a red giant branch star which had a high mass loss rate (e.g., Kilic et al. 2007). In addition, binary stars could also appear to be white dwarfs of masses < 0.45 M , for example unresolved white dwarf binaries such as double degenerate stars, look like one single overluminous star. With the same T eff , this means bigger radius, which means smaller mass according to the massradius relation of degenerate objects Bédard et al. 2017;Kilic et al. 2020). This phenomenon could mimic a white dwarf of mass < 0.45 M .
Higher-mass white dwarfs (> 1.1 M ) should also be viewed with caution. These stars are potentially magnetic white dwarfs and may have come from a merger (Ferrario et al. 2015), which would make the single star evolution assumption fail, and therefore the models would not reproduce the true evolution of the object.
The typical structure of a white dwarf consists of a C/O core with a thin shell of He, surrounded by a thin layer of H. Both the core and the surface layers have an important role in how the white dwarfs cool: the thermal energy is stored in the core while the heat outflow is logg/logg = 1% logg/logg = 5% Figure 10. Median relative uncertainty in the total age calculated using the Bayesian method in wdwarfdate. The median is calculated over all T eff and log g results with the same input uncertainties, and the uncertainty in the median is calculated as the standard deviation. In the left panel, the relative uncertainty in T eff is fixed at 1% and 10%. In the right panel, the relative uncertainty in log g is fixed at 1% and 5%. . Limits of T eff and log g for wdwarfdate: In brown are the values for which the code does not estimate any parameters because they are outside of the allowed range by the cooling tracks. In green-blue five-point stars are the values for which only a final mass and cooling age are estimated using the fast-test method, but not the rest of the parameters because of the limits where the MIST-based Cummings et al. (2018)  regulated by the opaque envelope (Hansen et al. 2004). For example, DA white dwarfs with thick H layers cool slower than non-DA objects because the H envelope is a better insulator (Fontaine et al. 2001). In addition, the chemical composition of the core will affect the heat capacity of the object, and therefore the cooling rate (Hansen et al. 2004). For instance, the estimated cool-ing age of a white dwarf assuming a pure C core is an upper limit, while assuming a pure O core is a lower limit, with differences of a few Gyr for the oldest objects (Fontaine et al. 2001). Although the core composition is still a source of uncertainty for white dwarf cosmochronology (Simon et al. 2015), it is necessary to assume one to estimate the cooling age and mass of a white dwarf with evolutionary models. wdwardate uses cooling tracks from Bédard et al. (2020) which assume a C/O core with the same proportions of C and O with the option of choosing between a thin or thick outermost H layer. The parameter estimation using wdwarfdate depends also on the relation between the progenitor mass and the white dwarf mass, meaning the IFMR. This is an empirical relation, which is improving as more precise measurements of the parameters of white dwarfs and their progenitors become available. The IFMRs included in wdwarfdate are similar to each other (See Figure 3), and we have seen no significant difference in the median values of the ages and masses estimated in this study when changing the IFMR. However, the IFMR should be chosen carefully when using wdwarfdate because the shape of the posterior PDF can be affected by it. Our recommendation among the available IFMRs is to use the MIST-based Cummings et al. (2018) relation, because this is the most complete and simplest relation to use with our method. Although Marigo et al. (2020) is more recent than Cummings et al. (2018), the difference is small between the relations, and the kink in the former complicates the sampling of the posterior without significant difference in the median mass or age of the white dwarf.
The main sequence age of the progenitor is obtained from a stellar evolution model, based on its mass, abundances and rotation, which can modify the estimation of the total age (Moss et al. 2022). We tested the parameter estimation for 50 randomly selected white dwarfs with 1500 K < T eff < 100, 000 K and 7 < log g < 9.3, with uncertainties of 10% and 1% respectively, using all the possible combinations of [Fe/H] = {−4, −1, 0.5} (dex) and v/vcrit = {0.0, 0.4}. For all the runs we used the MIST-based Cummings et al. (2018) IFMR, and the cooling sequences for DA white dwarfs. We compared the results to the parameters estimated using [Fe/H] = 0, and v/vcrit = 0.0 to estimate the effect of using different stellar evolution models. We found no significant change in the estimated values of the final mass or cooling age, which is expected because we did not modify the parameters of the white dwarf. For the initial mass we found a small difference of up to 0.16 M when [Fe/H] = −4 for both cases of v/vcrit. For the main sequence age and total age we found differences of up to 0.7 Gyr when [Fe/H] = −4, again not changing significantly with v/vcrit. This shows that the parameters of the progenitor star need to be chosen carefully given that they can affect the results. In case there is no information of the parameters of the progenitor star, it is reasonable to choose solar metallicity ([Fe/H] = 0), taking into account the possible variations described above.
The last aspect to take into account when running wdwarfdate is the choice of priors. Currently the code assumes a constant SFH and the IMF to be m −2.3 i when using the Bayesian mode (see Section 2.3), and it does not allow changes in these priors. As discussed in Sections 2.4 and 3.3, when the uncertainty in the input parameters is large, the likelihood will be less constrained and the priors will have a bigger effect, making the initial and the final masses smaller, and therefore affecting the estimated ages as well. In this case, the user could consider running the fast-test method instead of the Bayesian method. In this alternative method, wdwarfdate will perform a Monte Carlo uncertainty propagation to estimate the white dwarf parameters (as it was described in Section 2.4) and will not be affected by the priors.

SUMMARY
In this work we presented wdwarfdate, a publicly available Python package which derives Bayesian total ages of white dwarfs from an effective temperature (T eff ) and a surface gravity (log g), assuming single star evolution. wdwarfdate obtains the probability distributions by explicitly integrating over the likelihood and priors to perform the normalization for the following parameters: mass of the progenitor star, cooling age of the white dwarf and ∆ m , which models the scatter in the IFMR. From these parameters, the code obtains probability distributions for the mass of the white dwarf, the main sequence age of the progenitor and the total age. The derived ages depend on the input parameters and on the models chosen when initializing wdwarfdate, which are summarize in Table 1 (See Appendix A for an example on how to use wdwarfdate).
We tested wdwarfdate on white dwarfs cluster members from Cummings et al. (2018) and Canton et al. (2021). In general we found good agreement between our results and theirs for the mass and age of the white dwarf as well as the properties of the progenitor star. In addition, we used wdwarfdate to estimate the ages of a set of 18 white dwarfs which are candidate co-movers with M dwarfs (Kiman et al. 2021), and therefore conform a group of low-mass star age-calibrators.
We also described a detailed analysis of the constraints on the code which need to be taken into account before using. In summary: • We found a typical uncertainty of 10% on the Bayesian total age for T eff and log g with uncertainties of 1%, and of 25% with input uncertainties of 10% and 1% for T eff and log g, respectively.
• Combining the limitation of the cooling tracks with the restrictions of the IFMR, the values for which wdwarfdate can derive Bayesian total ages are 1, 500 K T eff 90, 000 K and 7.9 log g 9.3 (cm s −2 ).
• Given that the cooling tracks assume single star evolution and C/O core for the white dwarfs, the best range of final masses to use wdwarfdate is 0.45 − 1.1 M , because objects outside this range are not likely to have evolved as a single star.
• The IFMRs included in our code are similar to each other and the results are not modified significantly from using one over the other. However we recommend using the MIST-based Cummings et al. (2018) IFMR, given it is one of the most complete studies of white dwarfs in clusters currently available, and has a relatively simple functional form.
• The choice of [Fe/H] can modify the main sequence age and total age by up to 0.7 Gyr, but if there is no information about the progenitor star available, solar metallicity ([Fe/H] = 0) is a reasonable assumption.
• When the uncertainty of the input parameters is large, the prior in the Bayesian method will affect the estimation of the age given that the likelihood is not well constrained. The user can consider using the fast-test method in these cases.
Although Gaia has increased the number of known white dwarfs by a factor of 10, we still need more white dwarfs which belong to known clusters with spectra to improve the calibration of the IFMR. This relation is one of the main limitations when estimating the total age of a white dwarf because the complete functional form of the IFMR is still not known with great accuracy, especially in the low-mass regime.
In future work we plan to include the estimation of T eff and log g from photometry in wdwarfdate, and extra priors in case the final mass or total age of the white dwarf are known, to improve the estimation of the other parameters.
To learn about how to use wdwarfdate, check out the source code (Kiman 2022) 4 and the online documentation 5 .

B. RELATIVE UNCERTAINTY IN TOTAL AGE.
To study the relative uncertainty in the Bayesian total age estimated with wdwarfdate we run a grid of 100 white dwarfs with 1, 500 K < T eff < 100, 000 K and 7 < log g < 9.3, with different uncertainties. In this study we used the Cummings et al. (2018) MIST-based IFMR, with the DA white dwarf cooling sequences and [Fe/H] = 0 and v/vcrit = 0 stellar evolution models. In Figure 12 we show the relative uncertainty in the total age (t tot ) for four cases of input uncertainties: ∆T eff /T eff = 1% and ∆ log g/ log g = 1%, ∆T eff /T eff = 10% and ∆ log g/ log g = 1%, ∆T eff /T eff = 1% and ∆ log g/ log g = 5%, and ∆T eff /T eff = 10% and ∆ log g/ log g = 5%. As discussed in Section 4, the biggest effect on the uncertainty of t tot comes from increasing the uncertainty of log g. In addition, we show that higher T eff are the most affected by the increase in the uncertainty of the input parameters. We repeated this analysis for non-DA cooling models and found no significant difference. Figure 12. Relative uncertainty in the Bayesian total age estimated using wdwarfdate. We include four cases of different input uncertainties to show the variation as a function of T eff and log g.