Measuring the initial-final mass relation using wide double white dwarf binaries from Gaia DR3 Open Access

ABSTRACT

The initial-final mass relation (IFMR) maps the masses of main-sequence stars to their white dwarf descendants. The most common approach to measure the IFMR has been to use white dwarfs in clusters. However, it has been shown that wide double white dwarfs can also be used to measure the IFMR using a Bayesian approach. We have observed a large sample of 90 Gaia double white dwarfs using FORS2 on the VLT. Considering 52 DA + DA, DA + DC, and DC + DC pairs, we applied our extended Bayesian framework to probe the IFMR in exquisite detail. Our monotonic IFMR is well constrained by our observations for initial masses of 1–5 M_⊙, with the range of 1–4 M_⊙ mostly constrained to a precision of 0.03 M_⊙ or better. We add an important extension to the framework, using a Bayesian mixture-model to determine the IFMR robustly in the presence of systems departing from single star evolution. We find a large but uncertain outlier fraction of 59 ± 21 per cent, with outlier systems requiring an additional |$0.70_{-0.22}^{+0.40}$| Gyr uncertainty in their cooling age differences. However, we find that this fraction is dominated by a few systems with massive components near 0.9 M_⊙, where we are most sensitive to outliers, but are also able to establish four systems as merger candidates.

1 INTRODUCTION

After depleting the hydrogen fuel in their cores, almost all (≃97 per cent) main-sequences stars end their lives by transitioning to the giant phase and shedding their outer layers in an intense stellar wind. This mass-loss that takes place during stellar evolution stops short of the compact central core, which goes on to become a white dwarf. For most white dwarfs, the interior is composed of carbon and oxygen, enveloped by a layer of helium, and finally by a thin but opaque layer of hydrogen at the surface. The most massive white dwarfs, i.e. those with masses >1.1 M_⊙, will instead have interiors dominated by oxygen and neon (Camisassa et al. 2019). Main-sequence stars with masses above about 8 M_⊙ (Weidemann & Koester 1983), undergo further nuclear burning during the giant phases, leading to the formation of more exotic stellar remnants via core-collapse supernovae.

It is intuitive to consider that the mass of a main-sequence star ought to be correlated with the mass of the resulting white dwarf. More precisely the (final) white dwarf mass can be considered a function of the (initial) main-sequence star mass. This function is known as the initial-final mass relation (IFMR). The IFMR essentially encodes the amount of mass-loss occurring during stellar evolution for stars of different masses (Bloecker 1995), and so is an important component for understanding a wide range of astrophysical topics, from stellar evolution itself, to the chemical enrichment of the Galaxy. At a more basic level, the IFMR is also required simply to establish the progenitor masses for studies on individual white dwarfs.

Despite its universal importance, the precise functional form of the IFMR cannot be simply determined. From a computational perspective, determining the IFMR is limited by the predicted mass-loss occurring during the giant branches of stellar evolution, which is extremely sensitive to input physics, such convective and rotational mixing (Marigo & Girardi 2007; Cummings et al. 2019), and nuclear reaction rates (Fields et al. 2016).

From the observational side, the IFMR can be measured empirically,¹ taking advantage of white dwarf cooling as an accurate clock (Mestel 1952; Fontaine, Brassard & Bergeron 2001). Essentially, by accurately measuring the effective temperature (T_eff) and mass of a white dwarf, comparison with white dwarf cooling models yields the length of time since the white dwarf left the tip of the Asymptotic Giant Branch (AGB), i.e. the white dwarf cooling age. As long as an independent clock can be sourced to constrain the pre-white dwarf (pre-WD) lifetime, then the initial mass can be determined from stellar evolution models. Thus, for a sufficiently large sample of white dwarfs, the relation between initial mass, and final mass can be determined.

The most widely used observational method to determine the IFMR is using white dwarfs in clusters (Weidemann 1987; Ferrario et al. 2005; Salaris et al. 2009; Williams, Bolte & Koester 2009; Casewell et al. 2009; Dobbie et al. 2009), relying on the fact that the age of a cluster can be measured from the main-sequence turn-off. Therefore, all white dwarfs within a given cluster (which will have the same total age as the cluster itself) can have their pre-white dwarf lifetimes and hence their initial masses determined empirically. It then becomes a simple matter of fitting some function to the initial-mass/final-mass pairs. The most recent, and most extensive IFMR based on cluster white dwarfs was presented by Cummings et al. (2018), having examined 79 white dwarfs among 13 different clusters, and covering the entire range of initial masses from 0.85 to 7.5 M_⊙. Even so, such cluster studies are limited by the precision with which initial masses can be determined, and the low number of objects typically found per cluster.

However, cluster white dwarfs are not the only avenue for exploring the IFMR. As an alternative technique, El-Badry, Rix & Weisz (2018) used the population statistics of the Gaia DR2 white dwarf sample, to constrain the IFMR, finding reasonable agreement with results from cluster studies. Catalán et al. (2008) established a different technique to infer the IFMR, examining white dwarfs in common-proper motion pairs (i.e. wide binaries) with main-sequence stars. Both components of the system can be assumed to have formed at the same time, but evolving independently due to their wide separation. Specifically, when the more massive primary undergoes stellar evolution, no common envelope occurs, and therefore there is no mass transfer on to the secondary. As with clusters the final masses are determined from spectral modelling of the white dwarf, with initial masses estimated by age-dating the main-sequence companions using high-resolution spectroscopy. More recently, Barrientos & Chanamé (2021) applied a similar technique to 11 white dwarfs with turn-off/subgiant companions whose ages can be determined more precisely than for main-sequence companions. While accurate initial masses can be obtained this way, very few white dwarfs are found in binaries with companions at this specific evolutionary stage.

Also employing wide binaries, Andrews et al. (2015) prototyped a new technique for constraining the IFMR, instead using double white dwarfs (DWDs). Once again, both components of the binary are assumed to have evolved separately. However, unlike in the previous examples, one cannot directly infer the pre-WD lifetimes or total ages. Nevertheless, the cooling ages of the two white dwarfs still provide two independent clocks, where their difference (Δτ_WD) should have the same magnitude but opposite sign as the difference in pre-white dwarf lifetimes (Δτ_preWD), i.e. |$\mbox{$\Delta \tau _\mathrm{WD}$}= -\mbox{$\Delta \tau _\mathrm{preWD}$}$|⁠. This can be used to limit the combination of initial masses that can have produced that binary. By modelling DWDs within a Bayesian framework, Andrews et al. (2015) were able to construct a posterior distribution for the IFMR given a set of 19 DWD pairs. This work predates the release of Gaia DR2, with the authors acknowledging the impending deluge of Gaia data would open the opportunity for a much higher fidelity IFMR to be determined.

In this work, we adopt and extend the DWD Bayesian approach, first introduced by Andrews et al. (2015) for fitting the IFMR. In Section 2, we introduce our spectroscopically observed DWD sample. In Section 3, we discuss our approach to spectral modelling combined with precise Gaia photometry and astrometry. In Section 4, we explain the Bayesian framework we have adopted to model the IFMR, as well as our extensions to the framework first introduced by Andrews et al. (2015). We show and discuss our resulting IFMR fits in Section 5, testing a variety of assumptions and comparing with other established results from different methodologies. Finally, we give our conclusions and areas for future work in Section 6.

2 DOUBLE WHITE DWARF SAMPLE

2.1 Sample Selection

Our goal was to create as large a possible sample of DWDs. To achieve this, we selected double white dwarfs from Gaia DR2 (Lindegren et al. 2018) by combining the catalogue of El-Badry & Rix (2018) with our own selection. We selected all pairs from the WD catalogue of Gentile Fusillo et al. (2019) with parallaxes and proper motions consistent within 3σ. We limited our sample to those DWDs with accurate parallaxes (π/σ_π > 10) and projected separations of less than 10 000 au. We also required small errors on proper motion, |$\sqrt{\sigma ^2_{\rm pm\_ra} + \sigma ^2_{\rm pm\_dec}} \lt 6$| mas. Our chosen selection was designed to supplement El-Badry & Rix (2018) with DWDs at larger distances, and those that have poor-quality photometry in Gaia DR2.

Our selection includes 424 DWDs, whilst the catalogue of El-Badry & Rix (2018) contains 375. There is only moderate overlap between the two samples, with 214 DWDs being common to both: 161 objects are found in El-Badry & Rix (2018), but not in our sample. The majority of these (120) are close systems that pass the looser parallax and proper motion cuts of El-Badry & Rix (2018), that are designed to accommodate orbital motion of the DWD. A further 27 are explained by the relaxed separation constraint (50 000 versus 10 000 au) of El-Badry & Rix (2018), and 14 are missing from our selection for unknown reasons. Our selection includes 210 DWDs which are not reported by El-Badry & Rix (2018). The vast majority of these are either beyond 200 pc, or do not pass the cuts placed on the quality of photometry in El-Badry & Rix (2018). Of those 210, 23 systems are missing from El-Badry & Rix (2018) for unknown reasons.

Our final sample consists of 585 high-confidence DWDs, albeit with a complex set of selection criteria. With the improved astrometry of Gaia DR3 (Gaia Collaboration 2021), 99.8 per cent remain in our sample, and none of the rejected objects were observed.

2.2 Observations

From our sample of 585 DWDs, we observed 90 systems with the European Southern Observatory (ESO) Very Large Telescope (VLT) Focal Reducer and Low Dispersion Spectrograph (FORS2) under programmes 103.D-0718 and 109.213B. A summary journal of observations is provided in Table 1. Runs A, B, and C were observed under programme 103.D-0718 on the nights of 2019 June 1, 2019 July 2–3, 2019 September 24, respectively. Runs D and E were observed under programme 109.231B on the nights of 2022 June 2–3 and 2022 September 21, respectively.

For runs A, B, and C targets were selected based on brightness and visibility. Objects with existing Sloan Digital Sky Survey (SDSS) spectra were not observed and we preferentially chose systems with large cooling age and mass differences,²

based on the values provided in Gentile Fusillo et al. (2019). Prompted by the large fraction of magnetic white dwarfs amongst the systems observed in these runs (see Section 2.3 for details), for runs D and E we simply selected targets at random from the visible DWD systems without SDSS spectra. A Hertzsprung–Russel diagram of the observed systems is shown in Fig. 1, with Gaia astrometry/photometry and spectral classifications given in Table A1.

All observations were taken using the 1200B + 97 grating with a 0.7 arcsec slit, giving a wavelength range of 3660–5010 Å at a resolution of R = 2000. The slit was oriented to place both white dwarfs on the slit simultaneously. The data were reduced using the esoreflex FORS2 pipeline v5.6.2 and standard procedures. Spectra were extracted using optimal extraction (Horne 1986) and a modification of the pipeline was made to prevent the extraction window of faint targets being reduced if a brighter target is nearby on the CCD. Wavelength calibration was performed using a combination of arc lamp images and the location of bright sky lines.

2.3 Systems rejected from fitting

Out of our 90 observed systems, there were 18 systems we chose not to fit at all as we did not believe we could derive meaningful atmospheric parameters (red triangles in Fig. 1). For 12 of these, this was because at least one component exhibited a strong magnetic field. These systems are WD J0059 − 2417AB, WD J0220 − 1532AB, WD J0224 − 4611AB, WD J0344 + 1509AB, WD J0902 − 3540AB, WD J1159 − 4630AB, WD J1834 − 6108AB, WD J2018 + 2129AB, WD J2023 − 1446AB WD J2047 − 8206AB, WD J2304 − 0701AB, and WD J2353 − 3620AB. Note that while WD J2023 − 1446AB contains a strongly magnetic DA, it is primarily excluded from fitting due to its baffling A component (Section 2.4). Furthermore, four spectra exhibited an unambiguous helium-dominated atmosphere (either due to helium lines or pressure broadened metal lines). These systems are WD J0225 − 1756AB, WD J1813 + 0604AB, WD J1836 − 5114AB, and WD J2355 + 1708AB. Finally two systems, WD J1310 − 3930AB and WD J2115 + 2534AB, contained white dwarfs with blue DC spectra which could plausibly be explained by either strong magnetism or helium dominated atmospheres. In total, this left 72 systems used to perform spectroscopic fits.

2.4 Notes on individual systems

2.4.1 WD J1953 − 1019AB

This system is part of a resolved triple white dwarf, first reported by Perpinyà-Vallès et al. (2019). Our observations refer to the A and C components from the discovery paper. In hindsight, it would have been ideal to observe all three components of this system to provide further constraints on the IFMR. That said, Perpinyà-Vallès et al. (2019) found that the B and C components (the inner binary) have almost identical atmospheric parameters. The potential for this system to provide a strong constraint on the IFMR is therefore limited compared to an idealized case, where all three components have substantially different masses and cooling ages.

2.4.2 WD J2023 − 1446AB

While we are unable to use this system for fitting the IFMR, both components of this wide binary are noteworthy. The B component has a spectral type DAH. The H β absorption line shows Zeeman splitting indicative of a 14 MG magnetic field. The A component, however, is extremely unusual and while resembling a hot DQB, an assignment of DX (Sion et al. 1983) is more accurate: Its position in the Gaia HR diagram indicates a mass of ≃ 1.2 M_⊙ (Gentile Fusillo et al. 2021), and a plethora of lines are apparent in the spectrum (Fig. 2), none of which we have been conclusively able to assign. While a strong magnetic field with a non-hydrogen atmospheric composition provides a plausible explanation for these observations, if we allow for a large blueshift of 1050 km s⁻¹ (not accounting for the presumably substantial gravitational redshift), some of the strongest lines can be assigned to C ii, and a few to He i (Fig. 2), with no evidence of Zeeman splitting observed. Nevertheless, around one third of the spectral features remain unassigned, and some predicted spectral lines are absent. While the proposed blueshift could be explained as a wavelength calibration issue, we found no evidence for this, and the A component shows no similar shift, despite being observed simultaneously. Furthermore, we obtained a second set of observations of this system in our third FORS2 run (run C), finding a similar spectrum with features occurring at the same wavelengths, though many with different strengths.

Both components have similar proper motions of about 17 mas yr⁻¹, implying a v_⊥ of 15 km s⁻¹ for the system. Crucially, the difference in v_⊥ is <1.7 km s⁻¹ at the 99th percentile, which at an on-sky separation of 2000 au rules out the possibility of a chance alignment. Given the somewhat convincing velocity shift, it is tempting suggest WD J2023 − 1446A may itself be part of a close binary with a hidden companion, where both epochs were by chance observed at the same extremum in orbital velocity. While a much cooler white dwarf could remain hidden by the flux of WD J2023 − 1446A, the short orbital period (required for an orbital velocity >1050 km s⁻¹) would lead to orbital smearing of the spectra. An orbital period of a few 10 h could instead be achieved with a stellar-mass black hole companion, but such a scenario is obviously contrived. Light curves we have obtained of this object show no variability above the 1 per cent level. As it has no bearing on our investigation of the IFMR, we refrain from speculating further on the nature of this peculiar system.

2.4.3 WD J1336 − 1620AB

This system was first observed during observing run A, and subsequently re-observed during run D. While unintentional, it provided the opportunity to determine the level of systematic uncertainty in our spectroscopic fits. With essentially four observations (two per binary component), we were able to constrain the |$\mbox{$T_\mathrm{eff}$}$| relative precision to |$0.9_{-0.4}^{+1.0}$| per cent, and the log g precision to |$0.043_{-0.017}^{+0.045}$| dex. These precisions can be used as priors in Section 4.3.

2.4.4 WD J2018 + 2129AB

This system appears to be a rare wide binary where both components are magnetic. In the A component, Zeeman splitting of the Balmer lines is observed, indicative of a ≃ 1 MG field. At a glance, the B-component spectrum resembles a DC spectral type, though closer inspection reveals broad wavy absorption bands characteristic of a strongly magnetic atmosphere with a field strength in the 100s of MG (see WD J1159 − 4630B, WD J1834 − 6108A WD J2047 − 8206A, WD J2304 − 0701A, or WD J2353 − 3620B for more obvious examples). To our knowledge this is the only known wide DWD where both components exhibit magnetism.

3 SPECTRAL MODELLING

3.1 DA model grid

To obtain atmospheric parameters for our double white dwarf sample, we calculated a two-dimensional grid of pure hydrogen model spectra using the Koester white dwarf model atmosphere code (Koester 2010, 2013). Effective temperatures (T_eff) were calculated from 3500 to 10 000 K in steps of 250 K, then up to 20 000 K in steps of 500 K, and finally up to 35 000K in 1000 K steps. The grid was evaluated for surface gravities (log g) from 7.00 to 9.50 in 0.25 dex steps (cgs units). All models were calculated with a convective mixing length parameter ML2/α (Tassoul, Fontaine & Winget 1990) of 0.8. This choice of 0.8 allows easy application of 3D corrections to our spectroscopic parameters (see below).

Each model spectrum was convolved to an instrumental resolution of 1.8 Å, as appropriate for our FORS2 grating and slit-width combination. For each combination of T_eff and log g, we determined the stellar radius, R_WD, using the mass–radius relation from Bédard et al. (2020) appropriate for white dwarfs with thick hydrogen layers. Since the Koester model fluxes are in units of 4 × Eddington-flux, we scaled each model in the grid by a factor |$\pi (R_\mathrm{WD}/1\, \mathrm{kpc})^2$|⁠. As an additional step, we applied gravitational redshifts to each model in our grid, using same mass–radius relation as before. In principle, this allows our fits to be sensitive to the difference in gravitational redshift for the two components of a binary (assuming negligible orbital reflex motion), with the systemic radial velocity necessarily included as an additional free parameter.

To generate fluxes for arbitrary T_eff and log g and wavelength, we used trilinear interpolation of the logarithm of the model grid fluxes. The use of the log-fluxes is necessary as wavelengths bluer than the peak have a strongly non-linear dependence on T_eff.

Since these models are one-dimensional, they require correction to the corresponding 3D parameters for the most physically accurate results. We adopted the corrections given by Tremblay et al. (2013), applying their equations (9) and (10), appropriate for the ML2/α = 0.8 used in our model atmosphere calculations.

3.2 Spectrophotometric fitting

To fit each DWD, we used not only the FORS2 spectra, but also the Gaia DR3 photometry to constrain the atmospheric parameters. We first fit each binary using a least squares approach, following this up with a Markov chain Monte Carlo analysis for more informative parameter distributions.

Inspecting the FORS2 spectra, it was clear that the flux calibrations were generally poor, though in a consistent manner for each DWD, i.e. unphysical curvature in the spectral continuum of the A component would also be observed in the spectrum of the B component (as a result of both stars being observed on the slit simultaneously). We therefore took advantage of these shared systematic uncertainties by fitting both components of the binary simultaneously. Specifically, we modelled the wavelength dependent flux calibration correction as a fifth-order polynomial. While we found this approach to be generally successful, we found that the relative fluxes can differ by a constant factor, e.g. if one component was slightly off-centre on the slit compared to the other. We therefore introduced a further free-parameter to re-scale the fluxes of the B-component by a constant amount. A demonstration of our approach to flux calibrations is shown in Fig. 3 for the system WD J0023 + 0643AB.

Initially, we performed a least squares fit to the spectra and Gaia photometry³ to constrain the flux correcting polynomial, the flux scaling factor on the B-component, the T_eff and log g for each component, and the systemic radial velocity. Furthermore, we included the parallax as a free-parameter in order to marginalize over the uncertainty of the Gaia measurement. We then used the best-fitting parameters and covariance matrix to initialize a Markov chain Monte Carlo (MCMC) fit via the python package emcee (Foreman-Mackey et al. 2013).

As part of the MCMC, we introduced two new free parameters. In a few cases, where one white dwarf is significantly brighter than its companion, the bright component may have an extremely high S/N ratio (in order to get sufficient counts for the fainter component). This may cause that spectrum to dominate the fit entirely. To avoid this, we introduced a relative noise floor parameter, η to our fits. Simply put, we increased the flux uncertainties according to |$\sigma _i^{\prime 2} = \sigma _i^2 + (\eta \times f_i)^2$|⁠, where |$\sigma _i^{\prime }$| and σ_i are the modified and modified flux uncertainties, and the f_i are the spectroscopic fluxes at wavelength i. Compared with the usual |$-\tfrac{1}{2}\chi ^2$|⁠, the addition of this noise-floor requires a slightly modified log-likelihood for each spectrum,

where χ² has the usual meaning and implicitly depends on the |$\sigma _i^{\prime }$|⁠. A uniform prior for η > 0 was used on this parameter. Similarly, we folded in a 1 per cent additional flux to the Gaia fluxes according to the Gaia documentation on the uncertainty on the absolute calibration scales.

Our sample of DWDs have Gaia parallaxes locating them at distances up to 250 pc, with a median of 90 pc. Therefore, it is reasonable to assume that the Gaia fluxes for some of the more distant objects will be affected by interstellar reddening. Indeed in our initial attempts at fits (which did not include reddening), we found some cases where a good fit to the spectra resulted in models that were bluer than the Gaia photometry, affecting primarily the more distant objects such as WD J1336 − 1620AB. To account for this, we included interstellar reddening, specifically E(B − V), as a an additional free parameter. For the MCMC fits, we used the Jeffreys prior, P(E(B − V)) ∝ E(B − V)^−1/2, which is naturally weighted towards lower values on a linear scale, but is also a proper prior as long as an upper-bound is provided. For each DWD, we queried the Bayestar19 3D extinction maps (Green et al. 2019) to determine the 99th percentile of E(B − V) at the specified distance of the system. This was then used as upper-bound on our reddening prior. For systems outside of the Bayestar footprint, we adopted E(B − V) = 10⁻⁴ as the upper limit if located within 100 pc of the Sun, or 0.15 otherwise. For the most nearby systems, where the queried 99th percentile may return exactly zero, we also set the upper limit at 10⁻⁴ for numerical stability.

For the priors on the remaining free parameters we used normal distributions on both log gs of |$\mathcal {N}(\log g; 8.0, 0.25)$| (units of dex) to approximate the white dwarf mass distribution in the absence of any other strongly constraining data, though in practice, even in the worst case we found the data constrained the log g to within ±0.11 dex. For the radial velocity r, we again used a normal distribution |$\mathcal {N}(r; 0, 200)$| (units of km s⁻¹) specifically to avoid r becoming unbounded when fitting DC + DC pairs. Finally, for the two T_eff, we used flat priors on a logarithmic scale to reflect the greater abundance of cool white dwarfs in volume limited samples compared to the hottest objects.

3.3 Spectroscopic results

For the MCMC fits themselves, for each system, we used an ensemble of 200 walkers. We found that most of these had fully burnt-in within about 2000 steps. Therefore, we chose to run for 5000 steps, using the final 1000 steps (thinned to every 5 steps) to report the final results. The best parameters for all 72 fitted binaries are shown in Table A2, including derived parameters (masses, cooling ages, and difference in cooling ages). Spectra, photometry and best-fitting hydrogen atmosphere models for all systems are shown in Figs A2–A32. These figures include systems that we did not fit (Section 3.3.1), but which are shown without models.

The white dwarf mass distribution of all objects used to perform IFMR fitting is shown in Fig. 4 (details of additional excluded systems are given in Section 3.3.1). These white dwarfs have a median mass of 0.626 M_⊙, with the distribution appearing broadly similar to white dwarf mass distributions from other works (Genest-Beaulieu & Bergeron 2014; Rebassa-Mansergas et al. 2015; Hollands et al. 2018).

3.3.1 Systems excluded for fitting the IFMR

The vast majority of our fits are consistent with both spectroscopy and photometry (Figs A2–A25), therefore giving us confidence that they are suitable to fit the IFMR. Even so, there are some systems where, for one reason or another, we do not trust our atmospheric parameters, and so exclude them from the following parts of this work.

Four systems contain magnetic white dwarfs with field strengths of a few MG.⁴ Due to Zeeman splitting, the resulting fits had synthetic photometry that were entirely inconsistent with the Gaia data. These systems are WD J0002 + 0733AB, WD J0240 − 3248AB, WD J1314 + 1732AB, WD J1636 + 0927AB, and WD J2259 + 1404AB. Note that WD J1535 + 2125B exhibits a weak <1 MG magnetic field, which does not appear to have affected the quality of its spectroscopic fit. Therefore, we see no reason to reject this system.

Many systems contain one or more DC white dwarfs. Some of these remain ambiguous as to whether they could have hydrogen or helium-dominated atmospheres. However, a few systems show severe disagreement with their best-fitting hydrogen-dominated model spectra. The white dwarfs in question are WD J1124 − 1234A, WD J1350 − 5025A, WD J1729 + 2916B, WD J1929 − 3000B, and WD J2122 + 3005A. In most of these cases, the fit attempts to reduce the strength of the hydrogen lines by reducing the |$\mbox{$T_\mathrm{eff}$}$|⁠, resulting in extremely poor agreement with the photometry. Regardless, these systems have been fitted with an incorrect atmospheric model, and so must be rejected.

As stated above, some systems containing one or more DCs remain ambiguous. These are generally DCs with |$\mbox{$T_\mathrm{eff}$}\lt 5000$| K, which would not show strong hydrogen lines regardless of their atmospheric composition. Therefore, we rely on the photometry to assess the quality of the fits. Indeed, some DA + DC systems such as WD J0109 − 1042AB show excellent agreement with the photometry of the DC component, suggesting that object does indeed have a hydrogen dominated atmosphere. However, WD J0007 − 1605A, WD J1827 + 0403A, and WD J1859 − 5529B⁵ show disagreement suggesting they instead likely have helium-dominated atmospheres. We therefore exclude these three systems.

Our sample also contains DC + DC pairs introducing further ambiguity, as the atmospheric parameters for both components are essentially determined only by the photometry. All of these DC + DC systems, e.g. WD J0104 + 2120AB, show good agreement with their photometry. However, we found that several of these DC + DC pairs contained one or more white dwarfs with low derived masses implying they cannot have formed via standard evolution, as the universe is not old enough to have produced such low-mass objects assuming single star evolution for both components. Therefore, such systems must be rejected on the basis that they cannot constrain the IFMR. These systems are WD J1014 + 0305AB, WD J1211 − 4551AB, WD J1557 − 3832AB, WD J1804 − 6617AB, WD J1827 + 0403AB, WD J1929 − 4313AB, WD J2122 + 3005AB, WD J2230 − 7513AB, and WD J2248 − 5830AB.

After rejecting these 20 problematic systems (WD J1827 + 0403AB and WD J2122 + 3005AB are rejected on account of multiple reasons), our final sample we used to fit the IFMR contained 52 systems (mass distribution shown in Fig. 4).

3.3.2 Comparison of results

Since all objects in our sample were drawn from the DR2 catalogue of Gentile Fusillo et al. (2019) and subsequently are present in the EDR3 catalogue of Gentile Fusillo et al. (2021), we sought to compare our combined spectroscopic + photometric fits with the photometric results of Gentile Fusillo et al. (2021). In Fig. 5, we show the difference between these sets of results for both T_eff and log g, though specifically for the 52 binaries (104 white dwarfs) selected to constrain the IFMR. In both panels, the uncertainties in T_eff/log g differences are generally dominated by the contribution from the Gentile Fusillo et al. (2021) results, since we use the same photometric data for our fits, but have the benefit of spectroscopy.

In the left panel, a small but statistically significant linear trend is seen in the difference in results as a function of temperature. Specifically, we find a gradient of 0.051 ± 0.007, or equivalently 11.8 ± 1.5 per cent dex⁻¹. Across the range of data, this means we find 4 per cent lower T_eff at |$\log (\mbox{$T_\mathrm{eff}$}/\mathrm{K})=3.65$| (≃ 4500 K), consistent results at 9300 K, and 5 per cent higher T_eff at |$\log (\mbox{$T_\mathrm{eff}$}/\mathrm{K})=4.4$| (⁠|$\simeq 25\, 000$| K). The results for the surface gravities are much more consistent, with the best-fitting linear trend only varying by about 0.005 dex across the observed range of data (Fig. 5).

The T_eff differences, while generally small are still worth investigating as small changes in T_eff will result in larger differences in cooling ages. The source of this discrepancy could result from differences in the model atmospheres used between our work (Koester models) and the work of Gentile Fusillo et al. (2021) (Tremblay models), the fact that we used spectroscopy and photometry rather than just photometry, or fitting methodology. As a test we tried fitting the Gaia photometry alone. In this case, the linear trend in T_eff differences disappeared, suggesting the addition of spectroscopy is responsible. We note that despite finding no linear trends, our photometric fits were on average 3 per cent lower in T_eff compensated by results 0.04 dex lower in log g, indicating some remaining difference between atmosphere models or fitting methodology.

Differences in spectroscopic and photometric stellar parameters are well documented (Bergeron et al. 2019). Despite these differences, we note that given the size of the errorbars, we are still generally within 1–2σ of the results from Gentile Fusillo et al. (2021). Given that our spectroscopy provides significantly more data with which to work with, it is likely that our results are more accurate. Even so, in Section 4.2, we demonstrate the ability of our IFMR model to account for systematic uncertainty in T_eff and log g parameters, and so this is not expected to affect our results on fitting the IFMR.

4 MODELLING THE IFMR

In this section, we explain our Bayesian model for investigating the IFMR using wide DWDs. This is essentially an adaptation of the framework first introduced by Andrews et al. (2015), though presented in a mathematically simplified way. Furthermore, we have expanded the framework, accounting for correlated and underestimated T_eff and log g uncertainties, as well as allowing for the possibility that some fraction of systems are outliers that violate the assumption of coeval stellar evolution. Our implementation of this model will be made available at https://github.com/mahollands/IFMR_DWDs.

4.1 Bayesian framework

Given a set of observed wide DWDs with measured parameters, we ultimately wish to determine the distribution of IFMRs consistent with our data. In Bayesian notation, this means our desired posterior distribution can be written as

where the terms on the right are the prior probability distribution on the IFMR, and the likelihood of obtaining our data given a specific IFMR, respectively. Since each of the DWDs represent independent observations, the likelihood in equation (2) can be written as a product of likelihoods for each double white dwarf (DWD_k),

More explicitly, the relevant measured observables for each DWD are the two white dwarf (final) masses (M_f1 and M_f2) and the difference in their white dwarf cooling ages, Δτ_WD. Furthermore these parameters have measurement uncertainties, which we encapsulate in a covariance matrix, Σ. Written out explicitly, the individual DWD likelihood in equation (3) can be written as

where we have dropped the index, as the DWD in question is no longer explicitly from a larger set of DWDs. All that remains is to determine the functional form of this likelihood, though how to do so is not immediately obvious.

If instead we had prior knowledge of the initial masses, M_i1 and M_i2, things become much clearer. The predicted final masses could then be determined directly from the initial masses and IFMR. The difference in cooling ages could also be predicted from the main-sequence lifetimes corresponding to these initial masses, recalling that |$\mbox{$\Delta \tau _\mathrm{WD}$}= -\mbox{$\Delta \tau _\mathrm{preWD}$}$| for coeval systems. Mathematically, the likelihood for a DWD with a given IFMR and given initial masses can be written as

where X is defined as

In principle this likelihood could be used directly to determine the IFMR, by constructing a posterior distribution where the initial masses per system are also free parameters to be sampled. In practice, however, for large sets of DWDs such as ours, the posterior distribution becomes too highly dimensional to sample in finite time.

To ensure that the posterior distribution has as few free-parameters as possible (i.e. the parameters defining the IFMR), we can instead marginalize over the initial masses by integrating over them, i.e.

recovering the desired likelihood in equation (4). Conceptually, the process of integration takes into account all possible combinations of M_i1 and M_i2 leading to a likelihood where the distribution of data depends only on the IFMR. Even more rigorously, the integrand in equation (7) can be multiplied by a prior-distribution on M_i1 and M_i2 to weight the distribution by initial masses that are more common, i.e. an initial mass function (IMF) as a prior. In our implementation, we use a Salpeter IMF (Salpeter 1955) with exponent α = 2.3, i.e. the high-mass part of the Kroupa IMF (Kroupa 2001), and apply this to both components.

In practice, this integration step must be performed numerically, presenting its own set of challenges. The integral must be calculated many times (specifically the number of times the likelihood in equation (3) is evaluated multiplied by the number of DWDs in the sample), therefore computationally expensive integration techniques such as Gaussian quadrature are not appropriate here. Furthermore, for some specific DWD and IFMR sample, only a small region of the M_i1–M_i2 plane will contribute any significant probability density, making the choice of integration limits difficult (in order to avoid integrating over regions of near-zero probability density).

The joint-distribution of final masses, however, is known, and so for any given IFMR, the corresponding region in the M_i1–M_i2 plane can be determined, subject to the condition that the IFMR is monotonic (and thus invertible). Therefore by performing a change of variables, we can instead perform the integration in the M_f1–M_f2 plane, transforming the integral in equation (7) into

where the Jacobian terms are simply the gradient of the inverse IFMR evaluated at M_f1 and M_f2, and the M_i are determined by plugging the M_f into the inverse IFMR. Since the joint distribution of M_f1 and M_f2 has already been measured from the data, we can use Monte Carlo samples from this distribution, restricting ourselves to a small area in the M_f1–M_f2 plane (and the equivalent area in the M_i1–M_i2 plane) to perform the integral in equation (8).

We achieve this numerical approach using the technique of importance sampling. In importance sampling, one can evaluate the integral of a function f(x) (where f may be multivariate) by drawing N samples, x_i, from a distribution P(x). The integral may then be approximated as

Therefore importance sampling essentially turns the integral into a weighted mean, where – assuming that P(x) has been chosen to be similar in shape to f(x) – the density of samples is highest close to where f(x) is maximized, and few samples are placed at values of x that contribute little to the integral of f. Thus, the sample weights are simply 1/P(x_i), and where P is normalised over the space x.

In our case, we draw samples of M_f1 and M_f2 from a multivariate-normal distribution using the mean and covariance matrix derived from our spectroscopic fits, which can then be used to evaluate the double integral in equation (8). We demonstrate this approach in Fig. 6, where we use a 3-segment mock-IFMR to calculate the likelihood for WD J2007 − 3701AB in the M_f1–M_f2 plane. The red points represent 300 random draws from the measured values of M_f1 and M_f2, and their covariance matrix. The likelihood is slightly narrower than the distribution of Monte Carlo samples, as it is further constrained by the measurement of Δτ_WD. Although this means some of the 300 samples have very low likelihood, evaluating equation (7) directly would require sampling the entire M_i1–M_i2 plane, with almost all samples close to zero likelihood. We reiterate that integration over the M_f1–M_f2 plane via importance sampling, while drastically improving the computational efficiency of calculating the likelihood, explicitly depends on the IFMR model being monotonic – an assumption for which some recent works suggest may not be justified (see Sections 5.4 and 5.5).

With a numerical recipe for calculating the integral in equation (8), we have everything we need to calculate the total likelihood in equation (3). All that remains to determine the posterior distribution on the IFMR, P(IFMR|DWDs) in equation (2), is the prior distribution on the IFMR. In the following subsections, we introduce extensions to this framework. The priors on the associated hyperparameters are defined therein. Specific choices of priors on the IFMR and their effects on the results are given in Section 5.

4.2 Dealing with outliers

While some systems in our sample (e.g. WD J2223 + 2201AB, Table A2) clearly cannot obey an IFMR (i.e. two stars that formed together, evolving independently), other outliers may be more subtle. It is therefore insufficient to cherry-pick a-priori a subsample of our double white dwarfs in order to fit the IFMR. Instead, we have implemented a mixture model (Hogg, Bovy & Lang 2010) that naturally accounts for the presence of outliers.

We first consider the likelihood given for a double white dwarf where both components have evolved independently as single stars, i.e. equation (5), which we now refer to as L₀. For systems that have not followed this evolutionary path, e.g. due to mergers or additional mass-loss, we make a simple modification to L₀, instead assuming that the affected data have been drawn from a distribution with a wider variance than implied by the measured Δτ_WD. Therefore in the covariance matrix, Σ, we make the substitution |$\sigma _{\mbox{$\Delta \tau _\mathrm{WD}$}}^2 \mapsto \sigma _{\mbox{$\Delta \tau _\mathrm{WD}$}}^2 + \sigma _\mathrm{outlier}^2$|⁠, where σ_outlier is a free parameter representing the additional uncertainty in the cooling age difference. We refer to this altered uncertainty as L_outlier. Finally, these separate likelihoods are combined into a single likelihood, L_total, where

and where f_outlier is a free-parameter representing the fraction of systems in our sample that are outliers. This combined likelihood is used in all fits presented in later sections. In terms of priors, we used a uniform distribution on f_outlier in the range 0 to 1, and a uniform distribution on σ_outlier in the range 0 to 13.8 Gyr.

4.3 Dealing with underestimated uncertainties

Uncertainties for spectroscopic fitting of DA white dwarfs are generally assumed to be around a 1–2 per cent in T_eff and a few 0.01 dex in log g (Liebert, Bergeron & Holberg 2005; Voss et al. 2007). However, when fitting high-quality spectra such as ours (especially with the additional constraint of Gaia photometry), the formal uncertainties on the results are often found to be much smaller than the numbers above. Our analysis in Section 5 is no exception with instances of T_eff uncertainties below 10 K and log g errors of a few 0.001 dex (Table A2).

The effect of ignoring such small uncertainties is that the fitted IFMR may become overly constrained by a few systems. Furthermore, in practice we found that this led to some objects that would consistently yield a likelihood (equation 4) of zero which would propagate through to the total likelihood (equation 3). In the analysis of Andrews et al. (2015), the authors folded in T_eff relative uncertainties of 1.2 per cent, and log g uncertainties of 0.038 dex according to Liebert, Bergeron & Holberg (2005). While we could indeed adopt these same external uncertainties and apply them to our analysis, we instead chose to estimate them as part of our fits to the IFMR, including the relative T_eff uncertainty and log g uncertainty as free parameters to be determined.

We implemented this in our model by adding these uncertainties in quadrature into the T_eff–log g covariance matrices for each system, at each step in the fit. These were then converted to mass-Δτ_WD covariance matrices via the mass–radius relations and cooling age calculations of Bédard et al. (2020).

To treat these extra parameters as part of our Bayesian framework also requires considering their prior distributions. With no other prior information, an uninformative Jeffreys prior can be used, i.e. P(σ) ∝ 1/σ. Fortunately, we do have such extra information, as WD J1336 − 1620AB was observed twice, with somewhat different results on the T_eff and log g, exceeding the internal fitting errors.

For N observations of a quantity x (⁠|$\ln \mbox{$T_\mathrm{eff}$}$| or log g) for both white dwarfs in a binary, we consider the case where the values of x are drawn from normal-distributions with separate means (μ_A and μ_B), but a shared dispersion σ. Therefore, the likelihood of our data x is

However, since we only care about the value of σ we can marginalize over μ_xA and μ_xB, yielding

an inverse gamma-distribution in terms of σ² where

and where 〈x_A〉 and 〈x_B〉 are the sample means. Setting N = 2 and multiplying by the same Jeffreys prior as before transforms the likelihood into a posterior distribution on σ,

which in turn can be used as prior distribution on the uncertainty parameters when fitting the IFMR.

For our systematic relative |$\mbox{$T_\mathrm{eff}$}$| uncertainty (σ_T), and systematic log g uncertainty (σ_g), our measured values of WD J1336 − 1620AB led to S_T = 1.577 × 10⁻⁴ and S_g = 2.978 × 10⁻³. Consequently, we found |$\sigma _T=1.1_{-0.4}^{+1.1}$| per cent and |$\sigma _g = 0.046_{-0.018}^{+0.046}$| dex (median±16/84th percentiles), with 95 per cent highest density intervals of 0.3–4.0 per cent and 0.015–0.172 dex, respectively.

The median values for repeat observations of WD J1336 − 1620AB, are close to those reported by Liebert, Bergeron & Holberg (2005) (see above). Of course, systematic uncertainties in T_eff and log g can arise not only from the observations, but can also be stem from uncertainty in the models. These priors only estimate uncertainty from the former, but by allowing these additional variances to be free parameters, we allow additional sources of uncertainty to be accounted for in the posterior distribution, should the data support it.

4.4 IFMR functional form and parametrization

Thus far we have referred to the IFMR as some function, but without specifying its form, nor how it is parametrized. A natural choice of function is piece-wise linear, since it is easy to enforce monotonicity, as well as having a well-defined inverse, and inverse gradient.

For an IFMR comprised of N linear segments, this requires N + 1 free parameters, though the specific parametrization is implementation specific. For example the first segment could be specified simply by a gradient and intercept, with subsequent segments specified only by a gradient. Most published IFMRs are given in this format or similar (e.g. Williams, Bolte & Koester 2009; Salaris et al. 2009; Cummings et al. 2018), with each linear segment expressed in y = mx + c form. Andrews et al. (2015) used an alternative formulation in their analysis, specifying the angle of each segment makes in the M_i–M_f plane, with a final free-parameter specifying the perpendicular distance the extrapolated first segment makes with the origin. In our analysis, we simply chose to pre-specify a number of initial mass values, making the corresponding final-mass values the free parameters, with a similar approach also adopted by El-Badry, Rix & Weisz (2018). This choice makes it easy to check that the IFMR is monotonic, easy to invert the IFMR, and simple to calculate arbitrary values between break points (via linear interpolation).

Rather than restrict ourselves to a few segments, with the initial-mass values of the break points as free parameters, we instead chose to use a much finer grid of fixed initial masses. This allowed us to sample the IFMR across the entire initial-mass range, thus permitting us to detect subtle features in the IFMR without imposing any expectation on where they should occur. Furthermore, this would allow us to detect which regions of the IFMR are well constrained by our data, and which regions remain uncertain.

After experimentation, we settled on an initial mass grid at 12 fixed values of initial mass (11 segments), with grid points placed at 0.75–1.5 M_⊙ in 0.25 M_⊙ steps, from 1.5 to 4.0 M_⊙ in 0.50 M_⊙ steps, from 4.0 to 6.0 M_⊙ in 1.0 M_⊙ steps, and with a final step at 8.0 M_⊙. This grid is used throughout for the different fits described in the following subsections.

4.5 Pre-white dwarf lifetimes

The final ingredient required to constrain the IFMR is a relation between pre-white dwarf lifetimes and its dependence on initial mass. For this purpose, we used model grids from MIST (MESA Isochrones & Stellar Tracks; Dotter 2016; Choi et al. 2016), which themselves make use of MESA models (Modules for Experiments in Stellar Astrophysics; Paxton et al. 2011, 2013, 2015, 2018). Specifically, pre-white dwarf lifetimes were determined from the zero age main sequence to the beginning of the Thermally pulsing AGB (TP-AGB), with Solar metallicity (i.e. [Fe/H] = 0), and v/v_crit = 0.4.

Metallicity is known to affect both the pre-white dwarf lifetime and amount of mass of mass-loss occurring for a star of a given initial mass (Meng, Chen & Han 2008), with higher metallicities resulting in increased mass-loss (Romero, Campos & Kepler 2015). In cluster studies, the progenitor metallicity can be determined from the metallicity of other cluster members, and indeed can be found to be non-Solar. While we cannot directly probe the metallicity for the progenitors of our DWD sample, these objects are disc stars and therefore can be expected to have had values close to Solar (Andrews et al. 2015). While some of the systems in our sample must have total ages of many Gyr, Rebassa-Mansergas et al. (2021) recently investigated the age–metallicity relation for a large sample of white dwarf-main sequence binaries, finding no dependence of metallicity on their ages. Therefore, assuming Solar metallicity can be justified for the DWDs in our sample. Of course, stars in the disc show a spread in metallicity either side of Solar; however, this represents a source of systematic uncertainty which we expect to be absorbed by the parameters in our extended Bayesian model.

4.6 Summary

Here, we briefly summarize the differences between our framework and the original framework proposed by Andrews et al. (2015). First, equation (5) accounts for parameter covariances between M_f1, M_f2, and Δτ_WD. This is particularly important for our fits, since the two spectra of each DWD were fit simultaneously assuming shared systematics. Even so, for systems with independent spectroscopic observations (i.e. M_f1 and M_f2 are determined independently), Δτ_WD will necessarily covary with both masses, since cooling ages are sensitive to white dwarf masses.

Secondly, we make an important change to the importance sampling step. Andrews et al. (2015) chose to perform this integration over the space of pre-WD lifetimes (see their equation 13). In our opinion, this choice adds unnecessary computational effort, as the Jacobian will become the product of twice as many derivatives. By performing the integration over initial masses, as we have done in equation (7), we eliminate the need to calculate the gradient of the pre-WD lifetime function. Subjectively, we also find the choice to integrate over initial masses to be conceptually simpler, since one generally imagines pre-WD lifetime to be a function of the initial mass, rather than the other way around.

The other main differences to our model are the inclusions of outliers and unmodelled uncertainty in the T_eff and log g of our fits. These are already discussed in detail in the preceding subsections.

5 RESULTS AND DISCUSSION

5.1 A Monotonic IFMR

For our main investigation of the IFMR, we chose a simple set of constraints. Chiefly, we required the IFMR to be monotonic. This choice allowed us to apply the importance-sampling integral substitution in equation (8), which relies on the IFMR having both a well-defined inverse and inverse derivative.

In terms of other constraints, we only required the physically realistic restriction of fractional mass-loss between 0 and 1, and that total system ages were <13.8 Gyr (within the uncertainties). Specifically, we did not enforce the final mass free-parameters (on the IFMR) to be below the Chandrasekhar limit. This is because, first, we would like to infer from our results at what initial mass is the Chandrasekhar limit reached (if at all). Secondly, our highest initial-mass point, at 8 M_⊙, is somewhat arbitrary, but since we require the IFMR to be monotonic, limiting the final-mass points below 1.4 M_⊙, will therefore push all preceding final-mass points down to lower values in an artificial way. Nevertheless, we will investigate the effect of such a choice in Section 5.2.

We sampled the posterior distribution of our model described in Section 4, again using emcee (Foreman-Mackey et al. 2013) to perform an MCMC. Due to the large number of free parameters (12 for the IFMR, plus 4 hyper-parameters), we used an ensemble of 1000 walkers. For the importance sampling, we used 10 000 samples per evaluation of the likelihood. We found that the fit converged in |$\approx 10\, 000$| steps though we continued to 15 000 steps to ensure a large number of samples could be obtained post burn-in.

The resulting IFMR is shown in Fig. 7, along with the mass-loss fraction occurring as a function of initial mass. This is compared with the IFMR from Cummings et al. (2018) (a wider comparison against other published IFMRs is performed in Section 5.5). All results are given in the Fit 1 column of Table 2. A corner plot of converged parameters is shown in Fig. A1. In general, we find good agreement between 1 and 5 M_⊙, with extremely tight constraints for initial masses of 1–2 M_⊙ with M_f uncertainties of 0.01–0.02 M_⊙. The steepening gradient occurs in broadly the same region as that found by Cummings et al. (2018), though we found this section of the IFMR to be somewhat wider being located at 2.5–4.0 M_⊙ (compared with 2.85–3.60 M_⊙ for Cummings et al. 2018). From 4 to 5 M_⊙, we find our median result is flatter than Cummings et al. (2018), though still within our widening 1σ contour.

For initial masses above 6.0 M_⊙, we have virtually no constraint on the IFMR due to the absence of white dwarfs with final masses significantly above 1.0 M_⊙ in our fitted sample. Specifically in the mass-loss panel, the MCMC samples are seen to have filled the entire available parameter space, subject to the condition of bounded mass-loss and a monotonic IFMR. Similarly, for initial masses below 1.0 M_⊙ we have effectively no data covering this region, and so again the fit is very poor. This is most obvious in the mass-loss panel, demonstrating that a wide range of mass-loss fractions are covered by our 2σ contour at an initial mass of 0.75 M_⊙.

For the hyper-parameters, firstly, we found a particularly large (though poorly constrained) f_outlier of 0.59 ± 0.21, and with a corresponding σ_outlier of |$0.70_{-0.22}^{+0.40}$| Gyr. The interpretation of f_outlier is discussed in more detail in Section 5.7, but essentially is constrained only by a few systems with high white dwarf masses. For the systematic uncertainty parameters, we found |$\sigma _{\mbox{$T_\mathrm{eff}$}}=1.36_{-0.63}^{+2.09}$| per cent, and |$\sigma _{\log g}=0.049_{-0.006}^{+0.007}$| dex. For |$\sigma _{\mbox{$T_\mathrm{eff}$}}$|⁠, this remains similar to the prior distribution evaluated in Section 4.3, though with a slightly higher median value, and larger dispersion. For σ_log g on the other hand, the distribution was found to be almost Gaussian, with small uncertainty, indicating this parameter is particularly sensitive to our complete set of data. These values are in close agreement to the |$\sigma _{\mbox{$T_\mathrm{eff}$}}=1.4$| per cent and σ_log g = 0.042 dex found by Genest-Beaulieu & Bergeron (2019). Converting these uncertainties into masses and cooling ages, we found that typical white dwarf mass errors of 0.03 M_⊙, and cooling age relative uncertainties of 5–20 per cent. Importantly, for these four hyper-parameters, the corner-plot (Fig. A1) shows no obvious covariance with the IFMR parameters, meaning small changes to the hyper-parameters do not strongly affect the IFMR itself.

5.2 Chandrasekhar mass limit

In the previous section, it is clear that our IFMR is poorly constrained for M_i > 5 M_⊙, due to an absence of data with M_f > 1.0 M_⊙. As a result, for M_i in the range 5–8 M_⊙, our posterior distribution easily fills the entire parameter space between our priors of a monotonic IFMR prior and mass-loss fractions between 0 and 1. Naturally, this means that many of the samples exceed the Chandrasekhar limit to an extreme degree.

Our code allows the Chandrasekhar-limit to easily be imposed as a prior on the M_f values at each fixed M_i. We therefore repeated our fit to the IFMR with this extra constraint enforced. Because of the reduced phase-space our MCMC could explore, the model converged faster, only taking about 6000 steps, though we ran to 10 000 for the same reasons as before.

The resulting best fit is shown in Fig. 8, and tabulated in column Fit 2 in Table 2. The four hyperparameters are virtually unchanged, with differences only seen for the IFMR final mass parameters. At the low-mass end, where M_i < 2 M_⊙, essentially no difference is seen between this new IFMR and the one presented in Fig. 7 (Fit 1). Of course, at the high mass end, where the IFMR was previously unconstrained, the 1σ and 2σ contours appear relatively tight with |$M_f=1.101_{-0.102}^{+0.172}$| at M_i = 8 M_⊙. However, as predicted, this has also had an effect on the intermediate range. Because we have imposed a monotonic IFMR, forcing the high-mass end to be below 1.4 M_⊙, has by extension also pushed down the final mass values at intermediate initial masses by about 1σ. The result is that this IFMR now disagrees with the Cummings et al. (2018) IFMR by more than 2σ between initial masses of 3.5–5.5 M_⊙. We therefore recommend using the IFMR (Fit 1) over this one (Fit 2), and only for initial masses below 5 M_⊙.

5.3 Monotonic mass-loss

Inspecting Fig. 7, it is apparent that over the well constrained range (1–5 M_⊙), mass-loss also appears to be monotonic with respect to initial mass. This trend only breaks down at higher masses where we are unconstrained by data, and so our MCMC samples the entire parameter space. We therefore introduced the additional assumption that not only the IFMR is monotonic, but so also is mass-loss as a function of initial mass. To investigate the effect of imposing this restriction, we repeated the fit from Section 5.1, the same number of walkers (1000) and marginalization samples (10 000) were used running for 10 000 steps. The only change compared to Fit 1 was that we imposed a prior constraining the mass-loss fraction to also be monotonic.

The resulting IFMR is shown in Fig. 9, and with the results tabulated as Fit 3. Compared to our results with non-monotonic mass-loss (Fig. 7), the IFMR is pushed down substantially at the high-mass end (even more so than in Fit 2), with a final mass value of |$1.06_{-0.07}^{+0.11}$| M_⊙ at 8.0 M_⊙. Therefore, an IFMR with this assumption is unable to produce near Chandrasekhar-mass white dwarfs through single-star evolution. Furthermore, between 3.0 and 6.5 M_⊙, we find more than 2σ disagreement with the IFMR of Cummings et al. (2018), implying that the constraint of monotonic mass-loss is unlikely. Consequently, it is possible that the mass-loss fraction does in fact peak somewhere between 5 and 8 M_⊙. Evidently this result, and that of Fit 2, demonstrates that care must be taken when choosing priors, as these can have an unexpectedly large effect on the results.

At the low-mass end, between 0.75 and 1.0 M_⊙, this IFMR is more tightly constrained than in either Fit 1 or Fit 2, as IFMRs with high mass-loss at low initial-mass are ruled out, with the median line appearing to tend towards 0. Although, very poorly constrained by our data (only our priors), this result is somewhat realistic as the lowest mass stars – all of which are still on the main sequence given their long lifetimes – will evolve directly into white dwarfs without significant mass-loss once their hydrogen fuel is depleted, due to the reduced ability or inability to ignite helium burning. Therefore, users of our results may wish to use the parameters from Fit 3 when estimating the initial masses of low-mass white dwarfs, though these results are clearly unrealistic at higher masses, where parameters from Fit 1 should be used instead.

5.4 Non-monotonic IFMR

Although we restrict ourselves here to monotonic IFMRs only, we note that our code does permit non-monotonic IFMRs to be fitted. However, this requires changing the domain of integration from the final–final mass plane to the initial–initial mass plane, where it is not possible to restrict the area of integration to a region smaller than the whole plane up to 8 M_⊙ for both components. Therefore in order to maintain precision of the integrals, the number of integration samples should be increased by a factor of 100 or more.

Given these computational requirements, we consider such a task beyond the scope of this work at present. However, we note that recent work by Marigo et al. (2020, 2022) suggest that the IFMR may indeed be non-monotonic over a small mass range (1.8–2.2 M_⊙), and may be responsible for the formation of carbon stars. A more detailed comparison with this result is given in the following section.

5.5 Comparison with other IFMRs

In Sections 5.1–5.3, we presented our different IFMR fits while providing a comparison only to the IFMR of Cummings et al. (2018). Of course many other works in recent years have presented other IFMRs, which we also wish to compare to. In Fig. 10, we compare our best-fitting results to a variety of published IFMRs over the range 1–5 M_⊙, where our results are well constrained. This includes several cluster based IFMRs (Weidemann 2000; Kalirai et al. 2008; Williams, Bolte & Koester 2009; Cummings et al. 2018; Canton 2018), the Gaia population based IFMR of El-Badry, Rix & Weisz (2018), and the fiducial model of Andrews et al. (2015) using DWDs. As with our own IFMR, many of these works also make the assumption of a monotonic IFMR (Andrews et al. 2015; El-Badry, Rix & Weisz 2018), or are comprised of only a few linear segments and so do not result in negative gradients over the wide initial–mass range covered by each segment. We also show the non-monotonic IFMR of Marigo et al. (2020).

Over the range of 2.5–5.0 M_⊙, we observe agreement with all comparison IFMRs, within our 2σ shaded region. Above initial masses of 3.0 M_⊙, the comparison IFMRs show even better agreement, with the majority within the 1σ shaded region. At these higher masses, we generally have less data to constrain the IFMR, so it is perhaps no surprise that we find the greatest agreement where our constraints are weakest. Even so, at M_i = 4.0 M_⊙, our final masses are still constrained to within ±0.033 M_⊙.

The largest disagreement is observed between 1.0 and 2.5 M_⊙, depending on the IFMR. For the cluster results of Cummings et al. (2018) and Canton (2018), we continue to find reasonable agreement over this range, noting that the latter consists of a single linear relation. Compared to the other cluster IFMRs (Weidemann 2000; Kalirai et al. 2008; Williams, Bolte & Koester 2009), we find final masses that are about 0.05 M_⊙ higher. However, these are the oldest studies we show; the two more recent cluster analyses are consistent with our results, which we take that as a sign that our results are accurate in this range.

We also find a similar level of disagreement with the result of El-Badry, Rix & Weisz (2018), which was based on modelling the distribution of DA white dwarfs in the Gaia DR2 Hertzsprung–Russel diagram. At the low-mass end, their IFMR has similar level of uncertainty reported, and closely follows the results of Kalirai et al. (2008) over the whole initial mass range shown.

Somewhat surprisingly, at M_i = 2.0 M_⊙, our results show the largest disagreement with Andrews et al. (2015), despite using the same general methodology. Even at higher masses, our results remain about 1σ above those of Andrews et al. (2015). Of course, our results have benefited from a larger sample size, with parallaxes and precise photometry from Gaia DR3 also providing more accurate spectroscopic parameters. Since we have taken great care to account for sources of systematic uncertainties in both the spectroscopic fits and in the IFMR fitting, we expect that our IFMR ought to be more accurate, and with realistic uncertainties.

Finally, we observe clear disagreement with Marigo et al. (2020), where their IFMR peaks near 1.8 M_⊙. This non-monotonic IFMR was invoked following the discovery of cluster white dwarfs with low initial masses (1.6–2.1 M_⊙). In particular, a white dwarf in NGC 752 was found with a lower final mass, but higher initial mass than members of R-147 and NGC 7789, implying a downwards trend in the IFMR, before turning upwards again. Since our fits assume a monotonic IFMR as a prior, our model cannot reproduce this feature. However, if the IFMR kink was present in our sample, we would expect our resulting IFMR to have been forced upwards through the middle of the triangular feature, rather than across its base. This is especially true, given that the 0.6–0.75 M_⊙ range in final masses is covered by a large fraction of our DWD sample (Fig. 4). Given that the work by Marigo et al. (2020) required observations of particularly old clusters, we can only speculate that those cluster members may not be representative of the disk stars in our analysis. As discussed in Section 5.4, in principle our code has the ability to fit non-monotonic IFMRs, but at greatly increased computational cost, opening up the possibility to investigate the reported IFMR kink with a larger number of DWDs in the future.

5.6 Determining progenitor masses of white dwarfs

We consider the common case of deriving progenitor masses from the masses of isolated white dwarfs and an IFMR. The simplest approach is to draw normally distributed samples of the white dwarf mass from a measurement and its uncertainty, and then use the inverse-IFMR to determine the corresponding initial mass. This can be repeated over all IFMR samples to account for uncertainty in the IFMR.

A better alternative is to consider the likelihood of a final mass, given its uncertainty, an initial mass, and an IFMR,

This likelihood can be turned into a posterior distribution on M_i by multiplying by a prior distribution on M_i. In the absence of data to constrain M_i, the only assumption that can be made is that the star was drawn from the mass distribution of all main sequence stars, i.e. the IMF. We therefore adopt the high-mass part of the IMF (valid for M_i > 0.5 M_⊙) as a prior, so that |$P(M_i) \propto M_i ^{-2.3}$| and

However, this distribution depends on a specific IFMR (as opposed to a distribution of IFMRs), and so we must marginalize over the entire space of IFMRs, i.e. samples from the posterior in Section 5.1. Numerically, this equates to taking the arithmetic mean over IFMR samples

where the summand must be normalized over the range of 0.75–8 M_⊙, to provide correct weighting over each |$P(M_i | \mathrm{IFMR}_k, M_f, \sigma _{M_f})$|⁠. Since the resulting distribution is one-dimensional, it is simple to draw M_i samples from the posterior without resorting to more complex methods such as MCMC.

Posterior distributions from our Bayesian approach are shown in Fig. 11, using a final mass uncertainty of 0.03 M_⊙ (typical for our DWD sample, when systematic uncertainties are accounted for). The distributions for M_f = 0.6 and 0.7 M_⊙ appear particularly well behaved, given that they correspond to the most tightly constrained part of our IFMR. While the distributions for M_f = 0.8 and 0.9 M_⊙ appear further from normal, they are still adequate to determine 16th, 50th, and 84th percentiles and draw conclusions about the initial mass. For M_f = 1.0 M_⊙, sharp features appear in the posterior, resulting from a lack of constraining data, and poor sampling by our IFMR break points. We therefore recommend using our IFMR results over the range of 0.53–0.95 M_⊙ (for lower masses, the implied median pre-WD lifetimes exceed the age of the universe).

In comparison to the simple approach described earlier, we found the Bayesian approach yields more robust results. For example, using M_f = 0.60 ± 0.03 M_⊙, the simple approach gives |$M_i=1.34_{-0.23}^{+0.67}$| M_⊙, which is similar to the |$1.36_{-0.24}^{+0.42}$| M_⊙ found for the Bayesian method (Fig. 11), but with slightly worse uncertainty. However, inspecting the distribution of M_i samples, the simple approach also shows spikes or discontinuities at each break point in the IFMR (similar to what is seen for M_f = 1.0 M_⊙ in Fig. 11), which are not present for the Bayesian approach in this mass range. We provide a tool within our python package for calculating these distributions from white dwarf masses.

We conclude by noting that equation (15) relies on the IFMR only in the forward direction, making no reference to the inverse-IFMR, and so this approach is equally valid for non-monotonic IFMRs. While the resulting M_i distribution will inevitably be multimodal, these modes will be weighted by the IMF.

5.7 Initial masses and outlier probabilities for double white dwarfs

Within our Bayesian framework for fitting the IFMR (Section 4), we marginalized over initial masses in order to keep the number of free parameters to an acceptable level. However, it is still possible to calculate these after the fit. However, we cannot simply apply the methodology from Section 5.6 to each component separately. Firstly, because M_i1 and M_i2 may have correlated uncertainties, secondly because we also have the constraint on cooling ages and their difference, and thirdly because we must consider the effect of the derived hyperparameters on our sample.

To account for these caveats, we modify equation (17) to

where θ_k is a vector of hyperparameters, and Σ is the covariance matrix of the DWD parameters. We use the same M_f priors and likelihood as in Sections 4 and 4.2. While the resulting distribution is now 2D, it is still possible to sample the M_f1–M_f2 plane on a fine grid and calculate the posterior probability at each point. M_i1–M_i2 samples can then be drawn from this grid weighted according to the posterior. We give results for the estimated M_i1, M_i2 values and their covariance in Table A3.

Similarly, we can determine the probability that a specific DWD from our sample is an outlier, P_outlier, after having fitted our IFMR. This requires converged samples of the IFMR and the hyper-parameters, θ, from our fit. For each sample, we calculate the outlier likelihood, L_outlier and the total likelihood L_total (defined according to equation 10), but where both likelihoods have been marginalized over the M_i1–M_i2 plane (using importance sampling from Section 4). Finally, to determine the probability that a specific system is an outlier, we simply require marginalizing over our samples of fitted parameters according to

where f_outlier is the hyper-parameter corresponding to the fraction of systems that are outliers. These P_outlier will have values close to 1, if the two components in a system are not coeval, i.e. |$\mbox{$\Delta \tau _\mathrm{WD}$}\ne -\mbox{$\Delta \tau _\mathrm{preWD}$}$| considering the uncertainties on both Δτ_WD and Δτ_preWD. This necessarily means that we require high precision on both Δτ_WD and Δτ_preWD in order to be confident that a system is an outlier, or that both its components are coeval. We provide our estimate for P_outlier in Table A3.

The first thing to notice is that the vast majority of systems have P_outlier within one per cent of the mean value of f_outlier. This is because for most systems, the uncertainty on the observed Δτ_WD dominates the 0.7 Gyr of additional variance found for outlier systems (particularly once the systematics on T_eff and log g are considered), and so our methodology cannot confidently conclude whether additional variance is required. This begs the question of how f_outlier should be interpreted, and why its value is so high.

We demonstrate where this outlier probability arises from in Fig. 12, where we show the outlier probabilities for each system as a function of the more massive white dwarf in the pair, M_f2, and the uncertainty on Δτ_WD with assumed T_eff and log g systematic uncertainties of 1.4 per cent and 0.052 dex included, respectively. Across most of the parameter space, the P_outlier are almost identical to the central value found for f_outlier (0.59). Only six systems stand out from the crowd. These are WD J0101 − 1629AB, WD J1215 + 0948AB WD J2115 − 0741AB, and WD J2223 + 2201AB with high P_outlier; and WD J1313 + 2030AB and WD J1336 − 1620AB with low P_outlier. Inspecting the parameters of the high outlier probability systems it is easy to see how this classification was made, for example WD J0101 − 1629AB, the more massive component has been cooling for about 490 Myr, whereas the lower mass companion has a cooling age of 720 Myr – the inverse of what is expected for two white dwarfs born at the same time, but evolving independently as single stars.

All six of these systems (confident outliers, and confident non-outliers) have one thing in common: the masses of the heavier component, which are all between 0.87 and 0.91 M_⊙. While objects resulting from mergers are expected in this mass range, this does not explain why the most confident non-outliers are also found here. It is certainly true that for the more massive systems, the Δτ_WD uncertainties appear substantially smaller, with 240 Myr for the worst of the six (WD J1215 + 0948AB). Instead we speculate that only for higher masses⁶ do we have enough precision on both Δτ_WD andΔτ_preWD to make any distinction between outliers and non-outliers. Therefore, the outlier fraction for all systems, f_outlier is largely down to this small subset of six leading to such a high (and poorly constrained) value of |$0.59_{-0.21}^{+0.22}$|⁠. Inspecting the bottom panel of Fig. 12, we note that a few systems also show slightly more distinct P_outlier (both above and below 0.59) for white dwarf masses near 0.7 M_⊙. At this mass range though, the balance is slightly towards systems with lower outlier probabilities.

For the four high P_outlier systems, even if there is a selection bias towards high outlier/non-outlier confidence near 0.9 M_⊙, we conclude that the massive components in those four systems are candidate merger products from within former hierarchical triple systems. To rule out an alternative hypothesis of binaries formed via captures in multistar interactions, we looked at the distribution of P_outlier against projected sky separation. All four high probability outliers are separated by only 400–1200 au making captured binaries extremely unlikely, and lending further weight to our explanation of mergers. The widest system in our fitted sample, WD J2351 − 1601AB at |$27\, 600\pm 700$| au, does have a slightly high P_outlier of 0.68, raising the prospect that this system could have formed via a capture event.

This approach to outlier probabilities provided one final opportunity to investigate the possibility of a non-monotonic IFMR as identified by Marigo et al. (2020). We generated synthetic binaries with initial masses and system ages covering a wide range of values. These were converted to final masses using the IFMR of Marigo et al. (2020), and the pre-WD lifetimes (and hence cooling ages) calculated from the same MIST models as before. We then determined the outlier probabilities using the same method as above (i.e. using our own IFMR results to make this determination). We found that the outlier probability could be maximized by placing the component with lower initial mass at the peak of the Marigo et al. (2020) IFMR, i.e. at |$M_i=1.85\, \mbox{$\mathrm{M}_\odot $}$|⁠, while also requiring short cooling ages to reduce the uncertainties on Δτ_WD. With the initial mass of the other component |$\ge 2.10\, \mbox{$\mathrm{M}_\odot $}$|⁠, we were able to produce outlier probabilities of ≃ 0.90 (though we were unable to produce the ≥0.97 probabilities found for our most extreme outliers), demonstrating that our methodology can be sensitive to outliers produced by a non-monotonic IFMR.

As a final exercise we selected systems from our sample with both component masses between 0.6 and 0.75 M_⊙ (Table A2), i.e. covering initial masses of 1.5–3.0 M_⊙ from the Marigo et al. (2020) IFMR, and with at least one of the cooling ages below 300 Myr. This selection yielded 12 systems, all with outlier probabilities (Table A3) close to the median of 0.59, with the highest two values at 0.62 (WD J0120 − 1622AB and WD J0309 + 1505AB). Since these systems do not have outlier probabilities well above the 0.59 median or approaching the 0.90 found from our synthetic systems, we conclude that our IFMR and DWD sample do not provide evidence for a non-monotonic IFMR.

6 CONCLUSIONS AND FUTURE WORK

We observed 90 DWDs using FORS2 spectroscopy. Of these, we were able to use 52 DA + DA, DA + DC, and DC + DC pairs to constrain the IFMR assuming a monotonic piecewise-linear functional form. Furthermore, we were able to do this including the contribution of un-modelled uncertainties in their T_eff and log g, while also establishing that several systems in our sample are outliers, potentially suggesting the mergers taking place within former hierarchical triple systems.

While the vast majority of our binary sample can be fitted well with hydrogen dominated model spectra, some of the white dwarfs clearly have helium dominated atmospheres, with the most obvious of these having DB and DZ spectral classifications. Additionally some DC white dwarfs, while at first ambiguous, turned out to incompatible with a hydrogen dominated atmosphere at the derived temperature (i.e. hydrogen Balmer lines would otherwise have been present). In the future, we intend to also include these systems by fitting models with helium dominated atmospheres instead. However, it should be noted that a first investigation of the IFMR for hydrogen-poor white dwarfs suggests that these objects have higher mass progenitors, i.e. that they follow a separate IFMR altogether (Barnett et al. 2021).

Looking further ahead, our sample is still relatively small, containing only 52 pairs used to constrain the IFMR, and only for final masses up to 1.0 M_⊙. This means that our IFMR is only well constrained for initial masses below 5–6 M_⊙. Above that range, our fits merely reflect our choice of priors. Of course, identifying ‘well behaved’ systems containing ultra-massive white dwarfs where both components are non-magnetic DAs, presents its own observational challenges. Fortunately, the entire Gaia white dwarf sample contains over 1200 wide pairs⁷ – more than an order of magnitude increase over the 90 we have observed so far. Over the next decade, many of these will be observed in multifibre spectroscopic surveys such as WEAVE, DESI, SDSS V, and 4MOST. Taking advantage of these upcoming systems will allow us to constrain the IFMR to even greater precision in the future.

A caveat to our implementation of the Bayesian framework is the time-complexity. Naïvely, equation (3) implies a linear time-complexity in the number of white dwarfs, N, i.e. |$\mathcal {O}(N)$|⁠. However, the importance sampling step to calculate the integral in equation (8) limits the relative precision of this quantity for each DWD. Therefore, to maintain a constant relative precision in the total likelihood (equation 3), the number of integration samples (per system) must also increase linearly with N, increasing the overall time complexity to |$\mathcal {O}(N^2)$|⁠. For our fitted sample of N = 52, this does not currently pose an issue, with the main IFMR fit taking a few days on a modern 10-core desktop machine, but will become problematic for samples with N > 100. Therefore, alternative integration techniques that can reduce this time complexity should be investigated in advance of the large number of systems that will be observed in the coming years.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their feedback which improved the quality of this manuscript. MAH acknowledges useful conversations with Pier-Emmanuel Tremblay regarding 3D corrections to 1D atmospheric models. MAH and SL were supported by grant ST/V000853/1 from the Science and Technology Facilities Council (STFC). SCP acknowledges the support of a STFC Ernest Rutherford Fellowship.

DATA AVAILABILITY

The spectroscopic data obtained in this work can be found from the ESO Science Archive Facility with Program IDs 0103.D-0718 and 109.231B. The MCMC chains from our IFMR fits will be made available alongside our python package, and also available upon reasonable request to the authors.

Footnotes

References

APPENDIX A: ADDITIONAL TABLES AND FIGURES