The origin of ultramassive white dwarfs: hints from Gaia EDR3 Free

ABSTRACT

1 INTRODUCTION

The kinematics of stars can encode important information about their history and origin. The dynamical history of the Galaxy is imprinted on the kinematics–age relations of stars in the solar neighbourhood, such as the trend for the velocity dispersion of stars from the local Galactic disc to increase with age due to disc heating, as seen in the age–velocity dispersion relation (AVR; Strömberg 1946; Roman 1950a, b; Wielen 1977; Nordström et al. 2004; Holmberg, Nordström & Andersen 2007, 2009; Seabroke & Gilmore 2007; Soubiran et al. 2008; Aumer & Binney 2009; Holmberg, Nordström & Andersen 2009; Casagrande et al. 2011; Sharma et al. 2014; Yu & Liu 2018; Raddi et al. 2022). Since the epoch of the HIPPARCOS mission (Perryman et al. 1997), the kinematics of stars in the solar neighbourhood have been extensively studied (e.g. Dehnen 1998; Dehnen & Binney 1998; Chereul, Crézé & Bienaymé 1999; Skuljan, Hearnshaw & Cottrell 1999; Nordström et al. 2004; Famaey et al. 2005; Antoja et al. 2008; Famaey, Siebert & Jorissen 2008), and the precise astrometry of the Gaia mission (Gaia Collaboration 2016), the successor to HIPPARCOS, has further advanced these studies (e.g. Bovy 2017; Kushniruk, Schirmer & Bensby 2017; Gaia Collaboration 2018b; Yu & Liu 2018; Rowell & Kilic 2019; Mikkola et al. 2022; Raddi et al. 2022).

Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2021) has been used to study the kinematics of white dwarfs in the solar neighbourhood, showing that the kinematic structure of local white dwarfs (e.g. Mikkola et al. 2022) is similar to that of local main-sequence stars (e.g. Antoja et al. 2012; Kushniruk et al. 2017; Gaia Collaboration 2018b). For example, Raddi et al. (2022) performed a detailed three-dimensional kinematic analysis of a sample of 3133 white dwarfs from Gaia EDR3 with radial velocity measurements from either Gaia or cross-matched spectroscopic observations, finding that their sample, which was mostly located within ∼300 pc of the Sun, consisted of ∼90–95 per cent thin disc stars, ∼5–10 per cent thick disc stars, and a few isolated white dwarfs and halo members. For the thin disc members, Raddi et al. (2022) determined their AVRs and found them to agree with previous AVR results that were determined from different samples of LAMOST-Gaia FGK-type stars without radial velocity information (Yu & Liu 2018).

Gaia has also revealed structure in the colour–magnitude diagram of white dwarfs in the solar neighbourhood, which has implications for our understanding of white dwarf cooling. Gaia Data Release 2 (DR2) revealed the so-called Q branch (Gaia Collaboration 2018a), an observed sequence of massive white dwarfs in the colour–magnitude diagram that is transversal to the theoretical cooling sequences for white dwarfs with mass ≳1 M_⊙ and coincides with the region of white dwarf crystallization (Tremblay et al. 2019). The phase transition from liquid to solid state during the process of crystallization results in the release of latent heat (van Horn 1968) and other energy associated with element sedimentation (Garcia-Berro et al. 1988; Segretain et al. 1994; Isern et al. 1997; Althaus et al. 2012) that slows the rate of white dwarf cooling, and this cooling delay results in a pile-up of white dwarfs that has been identified with the Q branch (Tremblay et al. 2019).

While a cooling delay due to white dwarf crystallization had been predicted for over 50 yr prior to Gaia DR2, Cheng, Cummings & Ménard (2019) presented evidence from Gaia DR2 for an extra cooling delay of ∼8 Gyr experienced by ∼6 per cent of ultramassive (1.08–1.23 M_⊙) white dwarfs on the Q branch that is not explained by standard crystallization models. By looking at the number density of white dwarfs as a function of the age inferred by white dwarf cooling models, Cheng et al. (2019) noted an excess number of ultramassive white dwarfs on the Q branch above the expected pile-up from the crystallization delay. Furthermore, Cheng et al. (2019) noted that for some ultramassive Gaia DR2 white dwarfs on the Q branch there was a discrepancy between the dynamic age, inferred from the transverse velocity using the AVRs of Holmberg et al. (2009) and Sharma et al. (2014) for the local thin and thick discs, and the photometric isochrone age, inferred from the photometry using white dwarf cooling models. Among ∼1.08–1.23 M_⊙ white dwarfs, they found a population moving anomalously fast relative to the velocity expected from the AVR for their photometric ages. Cheng et al. (2019) argued that this discrepancy could only be explained by an extra cooling delay, in addition to the delays due to crystallization and double white dwarf mergers, and suggested modifications to the treatment of  ²²Ne sedimentation in white dwarf cooling models as a possible mechanism. This prompted many efforts to explain this cooling anomaly through modifications to white dwarf cooling models (Bauer et al. 2020; Blouin et al. 2020; Caplan, Horowitz & Cumming 2020; Horowitz 2020; Blouin & Daligault 2021; Blouin, Daligault & Saumon 2021; Camisassa et al. 2021; Caplan et al. 2021).

In our recent work (Fleury, Caiazzo & Heyl 2022, hereafter Paper I), we used the improved Gaia EDR3 data to re-investigate this cooling anomaly by analysing the number density distribution of the photometric cooling ages for massive (0.95–1.25 M_⊙) white dwarfs identified in the white dwarf catalogue of Gentile Fusillo et al. (2021). We considered a variety of publicly available white dwarf cooling models, including both the oxygen–neon core models of Camisassa et al. (2019) and carbon–oxygen core models of Bédard et al. (2020) with different envelope thicknesses. These are all standard cooling models, by which we mean that they do not attempt to reproduce the anomalous ∼8 Gyr cooling delay proposed by Cheng et al. (2019). For each set of cooling models considered, the masses and ages of the white dwarfs were determined from the Gaia photometry, and the resultant cooling age distributions were compared to the expectation from the time-varying star formation rate observed by Mor et al. (2019) for Gaia DR2 main-sequence stars. Mor et al. (2019) discovered a star formation burst in the local Galactic disc 2–3 Gyr ago, which, if a uniform white dwarf birthrate was assumed, would appear as an excess number of massive white dwarfs around this age produced through single stellar evolution. In Paper I, it was indeed found that, under the assumption of a uniform birthrate, there was an apparent excess of white dwarfs both along the Q branch and below it, coinciding with the burst of star formation seen in the star formation history of main-sequence stars.

As part of the detailed analysis in Paper I, the sample was further subdivided into three equally spaced mass bins, and the number density distribution of photometric cooling ages according to each set of cooling models was statistically compared to the Mor et al. (2019) star formation rate. For the two lightest mass bins, 0.95–1.05 and 1.05–1.15 M_⊙, it was found that standard cooling models could produce cooling age distributions that were statistically consistent with the expected distribution from the star formation rate observed by Mor et al. (2019). For the most massive bin, 1.15–1.25 M_⊙, it was found that the cooling age distribution was well fitted by a linear combination of the distribution expected for single stellar evolution products (taken to be proportional to the star formation rate) and the distribution expected for the product of double white dwarf merger products [taken to be the convolution of the star formation rate and the merger delay time distribution calculated by Cheng et al. (2020) using population synthesis simulations] when approximately 40–50 per cent of the 1.15–1.25 M_⊙ white dwarfs that formed over the past 4 Gyr were produced through double white dwarf mergers. From this analysis, it was found in Paper I that the photometric cooling age distribution of ultramassive white dwarfs could be explained by accounting for the time-dependent star formation rate and the presence of a large fraction of merger products among 1.15–1.25 M_⊙ white dwarfs, without needing to invoke an anomalous ∼8 Gyr cooling delay.

In this work, we follow up the work of Paper I with a kinematic analysis of the transverse motions of the same sample of 0.95–1.25 M_⊙ white dwarfs in the solar neighbourhood that was considered in Paper I, subdivided into the same mass bins of 0.95–1.05, 1.05–1.15, and 1.15–1.25 M_⊙ and with masses and cooling ages determined using white dwarf cooling models found in Paper I to produce photometric cooling age distributions consistent with the expectation from the star formation history. For each of these mass bins, we estimate the transverse velocity dispersion as a function of age using different estimators and analyse the distributions of the separate transverse velocity components for several age ranges. We compare the empirical distributions to the expectation from the AVRs measured for main-sequence stars in the local thin disc and find a population of anomalously fast-moving ultramassive white dwarfs on the Q branch with kinematic features that are not explained by an extra cooling delay in white dwarfs originating from the local thin disc.

The structure of this paper is as follows. The data and models used in this work, along with the formal procedure for determining transverse motions, are described in Section 2. The results are presented in Section 3, and the implications for the origin of anomalously fast-moving ultramassive white dwarfs on the Q branch are discussed in Section 4. Finally, the results and discussion are summarized in Section 5.

2 METHODS

2.1 Data sample

The completeness of this 200 pc sample as a function of M_G was analysed in Paper I and found to be volume limited for white dwarfs with M_G ≥ 15. The limiting volume as a function of magnitude, V_lim(M_G), was determined in Paper I using a variant of the Schmidt (1968) estimator as follows. The cumulative number distribution as a function of sampling volume was constructed for magnitude-binned sub-samples, and V_lim was determined for each magnitude bin by finding the volume above which the distribution began to deviate from the linear relation expected for a complete sample (and which was realized at smaller volumes). The limiting volume for each white dwarf was then calculated by linearly interpolating between these magnitude-binned V_lim values as a function of M_G, taking the limiting volume for magnitudes M_G ≥ 15 to simply be the total volume, V_max, of the 200 pc sample (i.e. the volume of a sphere with a radius of 200 pc).

2.2 Cooling and atmosphere models

Using the publicly available wd_models package provided by Cheng (2020),² we use white dwarf cooling models to determine the photometric age and mass of each white dwarf from its absolute G-band magnitude, M_G, and colour, G_BP − G_RP. To determine masses and cooling ages from these photometry measurements, the atmosphere composition of the white dwarf is an important factor to consider. As described in detail in Paper I, we classified each source in our sample as having either a pure-H, pure-He, or mixed atmosphere based on the best-fitting results from Gentile Fusillo et al. (2021), as determined by the chi-squared values provided in the catalogue for their model fits: chisq_H, chisq_He, and chisq_mixed. Each source was classified as having an atmosphere composition corresponding to whichever model gave the smallest chi-squared value. Sources for which all three chi-squared values were empty were assumed to have pure-H atmospheres.

Each set of models from wd_models considered in Paper I included both H-atmosphere and He-atmosphere models. The evolutionary models were combined with the corresponding publicly available Montreal group synthetic colours³ for pure-H and pure-He atmosphere models with Gaia EDR3 band-pass filters to convert the photometry measurements of each source into mass and cooling age values. The synthetic colours map the magnitude and colour to the effective temperature and surface gravity, respectively, T_eff and log g, and the evolutionary models map T_eff and log g to cooling age and mass. Models of the appropriate atmosphere composition were used for each source as determined by the atmosphere classification from the Gentile Fusillo et al. (2021) catalogue described above. He-atmosphere models were used for sources classified as having pure-He or mixed atmospheres, while H-atmosphere models were used for sources classified as having pure-H or unknown atmospheres.

In our previous work (Paper I), we assessed the consistency of the photometric cooling age distribution determined using various models available through wd_models with the star formation history observed for Gaia DR2 main-sequence stars (Mor et al. 2019). In this work, we use only the thick Montreal white dwarf cooling models (Bédard et al. 2020). These models gave the best-fitting cooling age distribution for the 1.05–1.15 and 1.15–1.25 M_⊙ mass bins in Paper I. For the 0.95–1.05 M_⊙ mass bin, the La Plata models (Renedo et al. 2010; Camisassa et al. 2019, 2017) and thick Montreal models were both determined to produce cooling age distributions consistent with the star formation history and yielded similar p-values, so for consistency we use the thick Montreal models for all three mass bins. Using the La Plata models instead of the thick Montreal models for the lightest mass bin produces similar results.

Since the publication of Paper I, the La Plata group has produced cooling models for ultramassive white dwarfs (≳1.15 M_⊙) with carbon–oxygen cores (Camisassa et al. 2022). These models implement more realistic physics in the core of the white dwarfs than the Bédard et al. (2020) models, such as element diffusion and abundance ratios determined by progenitor evolution simulations, though there is still a great deal of uncertainty regarding the abundances expected for carbon–oxygen white dwarfs (e.g. Salaris et al. 2010; De Gerónimo et al. 2017; Wagstaff, Miller Bertolami & Weiss 2020; Althaus et al. 2021). While an analysis of the photometric cooling age distributions determined by the Camisassa et al. (2022) models, particularly in comparison to the star formation rate, as was done in Paper I is worthwhile, our choice of cooling model will not change the salient results presented in this work. Using the Camisassa et al. (2022) models instead of the Bédard et al. (2020) models would shift the particular ages at which the key features of the velocity dispersion appear, but it would not remove these features. Likewise, using oxygen–neon core models (e.g. Camisassa et al. 2019) instead of carbon–oxygen core models would shift the inferred cooling ages and, to a lesser extent, masses of our white dwarf sample, but it would not affect the values of the velocities and therefore would not eliminate the key features presented in this work.

2.3 Kinematics

2.3.1 Streaming motion

The observed velocity of a star in the solar neighbourhood consists of the local (average) streaming velocity, which includes both the average motion due to Galactic rotation and reflex solar motion relative to the local standard of rest, plus the peculiar velocity of the star due to random motion. In this work, we are interested in the random motion of white dwarfs about the local mean, so we correct the observed velocities for both Galactic rotation and solar motion.

We calculate the streaming velocities using values reported in the literature for both the solar motion and the Oort constants. For the solar motion, we use the peculiar solar velocity measurements of Wang et al. (2021), with values of U₀ = 11.69 km s⁻¹, V₀ = 10.16 km s⁻¹, and W₀ = 7.67 km s⁻¹. For the Galactic rotation, we assign the Oort constants the values measured by Bovy (2017) using Gaia DR1 main-sequence stars from the Tycho–Gaia Astrometric Solution catalogue, with values of A = 15.3 km s⁻¹ kpc⁻¹, B = −11.9 km s⁻¹ kpc⁻¹, C = −3.2 km s⁻¹ kpc⁻¹, and K = −3.3 km s⁻¹ kpc⁻¹.

2.3.2 Random transverse motion

Since most sources in Gaia EDR3 do not have radial velocity measurements, we perform our analysis of the kinematics using only the transverse motion. Following the procedure of Dehnen & Binney (1998), we relate the transverse motion to the full three-dimensional motion through use of a projection operator |$\mathbf{ \boldsymbol{ \mathsf{ P}}}$| that projects the three-dimensional space velocity |$\mathbf {\boldsymbol{ v}}$| onto the celestial sphere.

Throughout this work, we compare the results we find for the white dwarf transverse kinematics relations with the expectation assuming the dispersions of their full three-dimensional velocities follow the AVRs determined by Holmberg et al. (2009) for main-sequence stars from the local thin disc. While more recent determinations of these AVRs exist (e.g. Yu & Liu 2018), they are similar to the AVRs of Holmberg et al. (2009) and do not change our results. As evidence for a cooling anomaly among ultramassive white dwarfs was presented by Cheng et al. (2019) with reference to the thin disc AVRs of Holmberg et al. (2009), using the Holmberg et al. (2009) AVRs facilitates comparison of our results with earlier work on this topic.

The sample for which the Holmberg et al. (2009) AVRs were determined was volume-complete to a distance of ∼40 pc (Nordström et al. 2004; Holmberg et al. 2007, 2009), so both the Holmberg et al. (2009) sample and our 200 pc sample are local enough that bias against populations with large scale height is not a major concern in comparing our results to the Holmberg et al. (2009) AVRs. This was verified by comparing the Holmberg et al. (2009) AVRs to the Yu & Liu (2018) AVRs for a sample of thin disc (metal-rich) stars with |z| < 270 pc, which we found to be similar enough that the main features in the comparison of our results to the AVRs were unchanged. The similarity of the Yu & Liu (2018) AVRs to the Holmberg et al. (2009) AVRs demonstrates the negligible impact of scale height bias for the results presented in this work.

3 RESULTS

3.1 Summary statistics

Fig. 1 shows various probes of the AVR for Gaia EDR3 white dwarfs in three different mass bins, 0.95–1.05 M_⊙ (blue), 1.05–1.15 M_⊙ (orange), and 1.15–1.25 M_⊙ (green). The masses and cooling ages were calculated from the observed photometry using the single stellar evolution white dwarf cooling models of Bédard et al. (2020) which were determined in Paper I to produce cooling age distributions that are statistically consistent with the star formation rate observed for Gaia DR2 main-sequence stars (Mor et al. 2019).

Figure 1.

Upper panel: Fraction of fast-moving massive white dwarfs for three mass bins. The fast-moving white dwarfs are defined as those with transverse velocity greater than 60 km s⁻¹. Middle panel: Transverse velocity dispersion accounting for uncertainties in proper motion. Lower panel: Scaled median transverse velocity (robust estimator of velocity dispersion). In the middle and lower panels, the solid black line traces the results for main-sequence stars from Holmberg et al. (2009) with the same sky distribution as the white dwarfs, and the dashed black lines trace the results from main-sequence stars including a merger delay described by the merger delay distribution of Cheng et al. (2020). The main-sequence curve (solid black) is the expectation if the 1.15–1.25 M_⊙ white dwarfs are all produced by single progenitors, while the mergers curve (dashed black) is the expectation in the limit that 100 per cent of those white dwarfs are produced by double white dwarf mergers.

Open in new tabDownload slide

The upper panel of Fig. 1 shows the fraction of fast-moving white dwarfs with v_t > 60 km s⁻¹ as a function of photometric cooling age. The empirical distributions, shown as histograms, were constructed using cooling age bins of 0.5 Gyr width over the range 0–4 Gyr. The fraction of fast-movers in each age bin was calculated using weighted counts for that bin, where each white dwarf was assigned a weight according to the procedure of Paper I to correct for the reduced sampling volume that ensures completeness.

The fraction of fast-moving ultramassive white dwarfs in Gaia EDR3 as a function of photometric cooling age shows a trend with cooling age that is similar to what Cheng et al. (2019) found for Gaia DR2 white dwarfs, though we see a less pronounced peak for young 1.05–1.15 M_⊙ white dwarfs on the Q branch than what was found by Cheng et al. (2019) for this bin. For the 0.95–1.05 M_⊙ white dwarfs, the fraction of fast-movers tends to increase with cooling age, which is the expected behaviour if the white dwarf velocities follow an AVR similar to what has been observed for main-sequence (G and F dwarf) stars in the thin disc (Holmberg et al. 2009). Of particular note is the distribution for 1.15–1.25 M_⊙ white dwarfs, for which the fraction of young fast-movers in the age range 0.5–1.5 Gyr is larger than the fraction in either the younger or older age bins. This corresponds to an excess of young fast-movers on the Q branch relative to the expectation based on the AVR observed for stars in the thin disc (Holmberg et al. 2009). While we see a hint of a similar excess in the 1.5–2.0 Gyr bin of the 1.05–1.15 M_⊙ white dwarfs, we do not see a clear signal of an excess of fast-movers within the uncertainties. This difference between our distribution and that of Cheng et al. (2019) for the 1.05–1.15 M_⊙ white dwarfs could be due to the difference in colours and proper motions between the Gaia EDR3 and DR2 measurements. The updated colour measurements for very blue objects in particular affect the masses inferred for those objects.

The middle and lower panels of Fig. 1 show estimates of the transverse velocity dispersion as a function of age using different estimators of σ_t, one based on the sample mean in each bin (middle panel) and the other based on the sample median (lower panel), the latter of which is more robust to outliers. Random variation in the observed values of v_t will arise both from the intrinsic randomness in the motion of the white dwarfs about the local mean (i.e. the peculiar motion of the white dwarfs) that we are interested in and from random measurement error. We correct for the latter effect in both methods of estimating σ_t from the observations. The random error for each value of v_t was calculated by propagating the standard errors and correlations provided by Gaia for the proper motions and parallax from which the observed transverse velocity (without corrections for Galactic rotation or solar motion) was calculated. The error was propagated using the expression |$e_{v_\mathrm{ t},\mathrm{(obs)}}^2 = \mathbf{ \boldsymbol{ \mathsf{ J}}} \ \mathbf{ \boldsymbol{ \mathsf{ \Sigma}}} \ \mathbf{ \boldsymbol{ \mathsf{ J}}}^\mathrm{ T}$|⁠, where |$\mathbf{ \boldsymbol{ \mathsf{ \Sigma}}}$| is the covariance matrix, |$\mathbf{ \boldsymbol{ \mathsf{ J}}}$| is the Jacobian matrix, and |$e_{v_{\mathrm t},\mathrm{(obs)}}$| is the random error for the observed transverse velocity.

We similarly determined the expected relation for the Holmberg et al. (2009) AVRs if all of the white dwarfs are the product of double-white dwarf mergers, which is shown as the dashed black curve in the middle and lower panels of Fig. 1. This corresponds to the limit in which the merger fraction is 100 per cent for white dwarfs that formed over the last 4 Gyr. Instead of using the photometric cooling age inferred directly using the white cooling models applicable for single stellar evolution to calculate the predicted σ_U, σ_V, and σ_W values from the AVRs, the relation for merger products was calculated using the expected true age after accounting for a merger delay time, with the expectation value calculated using the merger delay distribution of Cheng et al. (2020). Delay time distributions associated with double white dwarf mergers were calculated by Cheng et al. (2020) using binary population synthesis simulations for resultant white dwarfs of various mass ranges; we use the distribution for 1.14–1.24 M_⊙ white dwarfs, which spans approximately the same mass range as our most massive bin and is the merger delay time distribution that was used in Paper I.

Double white dwarf mergers can produce photometrically young white dwarfs with transverse velocities faster than what would be expected for single stellar evolution products of the same photometric age because the true age of the white dwarf merger products (accounting for the lifetime of the merger progenitors) is typically older than the age inferred from photometry (Temmink et al. 2020). Estimates of the merger fraction from both simulations (e.g. Bogomazov & Tutukov 2009; Temmink et al. 2020) and data (e.g. Cheng et al. 2020; Kilic et al. 2021, 2023; Fleury et al. 2022) indicate that double white dwarf merger products constitute an appreciable fraction of ultramassive white dwarfs, with estimates ranging from ∼15 per cent to over 50 per cent. The population synthesis simulations of Temmink et al. (2020) indicate that 30–50 per cent of white dwarfs with mass ≥0.9 M_⊙ form through binary mergers, with about 45 per cent of these mergers occurring between two white dwarfs, while the population synthesis simulations of Bogomazov & Tutukov (2009) indicate that over 50 per cent of white dwarfs with mass ≥1.1 M_⊙ are produced through double white dwarf mergers.

Empirically, based on the magnetism, kinematics, rotation, and composition of the 25 most massive white dwarfs identified by Kilic et al. (2021) in the Montreal White Dwarf Database 100 pc sample, Kilic et al. (2023) found a merger fraction of |$56^{+9}_{-10}$| per cent among ultramassive white dwarfs with mass ≳1.3 M_⊙. Theoretical results indicate that nearly half of the identified merger products would have been produced through the double white dwarf merger channel (Temmink et al. 2020), and over half of the merger products identified in Kilic et al. (2021, 2023) showed kinematic and/or rotation features indicative of double white dwarf mergers in particular. Cheng et al. (2020) estimated a double white dwarf merger fraction of about 20 per cent for Gaia DR2 white dwarfs within 250 pc with mass in the range 0.8–1.3 M_⊙ and a larger fraction of about 35 per cent for the subset in the mass range 1.14–1.24 M_⊙ (assuming carbon–oxygen cores to determine the mass). For our Gaia EDR3 white dwarf sample within 200 pc in the 1.15–1.25 M_⊙ mass bin, it was found in Paper I that the photometric cooling age distribution was well fitted when the double white dwarf merger fraction was about 40–50 per cent.

The expectation for the limiting case of 100 per cent double white dwarf merger products is shown in Fig. 1 as a simple check of the feasibility of a double white dwarf merger delay explanation of the fast-movers. The lower panel of Fig. 1 shows that a sufficiently large fraction of merger products among 1.15–1.25 M_⊙ white dwarfs can reach the excess observed for the median estimate of the transverse velocity dispersion, which is representative of the (more slowly moving) bulk of the white dwarfs as the median estimator is robust to the outlying fast-movers. However, even in the extremum limit that 100 per cent of the 1.15–1.25 M_⊙ white dwarfs are the product of double white dwarf mergers, the middle panel of Fig. 1 shows that the excess observed for the mean estimate of the transverse velocity dispersion, which is sensitive to the outlying fast-movers, cannot be achieved. In keeping with Cheng et al. (2019), we thus find that double white dwarf mergers (in the local Galactic disc) cannot be the sole explanation of the fast-movers.

3.2 Three-dimensional velocity distributions

The behaviour of the AVRs that we measure in Section 3.1 for the transverse velocities aligns with the behaviour of the observed fraction of fast-movers as a function of cooling age. For both the non-robust and robust estimators of transverse velocity dispersion, we see that the lightest (0.95–1.05 M_⊙) white dwarfs tend to become more dispersed with age, as inferred from the fraction of fast-movers and expected from the AVRs of main-sequence stars. The 1.05–1.15 M_⊙ white dwarfs mostly follow this trend as well, but with a small peak for cooling ages of ∼1.0–2.0 Gyr. The most massive (1.15–1.25 M_⊙) white dwarfs show a prominent peak for cooling ages of ∼0.5–2.0 Gyr. The peak in transverse velocity dispersion of ultramassive white dwarfs inferred using the non-robust estimator (mean) is much larger than that inferred using the robust estimator (median), suggesting that the fast-movers in these age bins are outliers or comprise a separate population that is much more dispersed than the main population.

To better understand the kinematic features of the population of young fast-movers, we consider the distribution of each Cartesian component of the transverse velocities in each mass and age bin. Fig. 2 shows quantile–quantile (Q–Q) plots comparing the sample distribution of transverse Cartesian velocity components U_⊥ (upper row), V_⊥ (middle row), and W_⊥ (bottom row) to the quantile function for a normal distribution with zero mean. The quantile function of a distribution is the inverse of the cumulative distribution; for an input percentile p, it gives the value of the random variable for which there is p probability of the random variable being less than or equal to the output value. For a 1D normal distribution with mean μ and variance σ², the quantile function is |$\mu + \sigma \ \sqrt{2} \ \mathrm{erfc}^{-1}\left(2 \left[1 - p \right] \right)$|⁠, where p ∈ [0, 1] is the percentile (i.e. cumulative probability) and erfc⁻¹ is the inverse complimentary error function. In Fig. 2, the measured value of the velocity component is used as the x-axis variable, while the y-axis variable is the value of the quantile function in units of standard deviation calculated from the weighted cumulative fraction. Written explicitly, the y-axis variable is |$y_i = \sqrt{2} \ \mathrm{erfc}^{-1}\left(2 \left[1 - p_i \right] \right)$| for the ith ordered data point (in order of increasing x) with p_i = ∑_{j ≤ i}w_j where j runs over all points for which x_j ≤ x_i. For reference, we also show quantile functions for normal distributions with zero mean and σ = 10, 20, 30, and 40 km s⁻¹ in order of decreasing slope as black lines. If the transverse velocity in a particular direction is normally distributed for a particular sample, the corresponding Q–Q plot in Fig. 2 should appear as a straight line. A non-zero mean will simply cause a shift of the distribution along the horizontal axis, as can be seen in the plots for V_⊥ due to asymmetric drift.

Figure 2.

Quantile–quantile plots of transverse velocity components. The vertical axis is the quantile function in units of standard deviation σ for a normal distribution with zero mean, calculated from the empirical cumulative fraction corresponding to the velocity component shown on the horizontal axis. The black lines correspond to normal distributions with σ = 10, 20, 30, and 40 km/s in order of decreasing slope.

Open in new tabDownload slide

Observations of main-sequence stars (e.g. Holmberg et al. 2009) indicate that the velocity dispersion of stars in the local thin disc increases with age. If the white dwarf velocities follow the AVRs of main-sequence stars (Holmberg et al. 2009), the slope of the distribution in Fig. 2 should get progressively less steep with increasing cooling age for a given mass bin and velocity component. For the age bins considered in Fig. 2, this effect should be more pronounced for younger ages due to the power-law form of the AVRs with exponents <1. This is the approximate behaviour seen for the lightest two bins, 0.95–1.05 M_⊙ and 1.05–1.15 M_⊙ (left and middle columns of Fig. 2). For the ultramassive 1.15–1.25 M_⊙ white dwarfs, however, we see that the curves for the 0.5–1.5 Gyr white dwarfs have a shallower slope and thus larger dispersion than what would be predicted by the Holmberg et al. (2009) AVRs, particularly for negative U_⊥ and V_⊥. The dispersions for these 0.5–1.5 Gyr white dwarfs are notably larger than the dispersions for the older 1.5–2.5 and 2.5–3.5 Gyr white dwarfs with mass in the range 1.15–1.25 M_⊙. Furthermore, for the 0.5–1.5 Gyr age bin of 1.15–1.25 M_⊙ white dwarfs, the white dwarfs with negative V_⊥ are much more dispersed in V_⊥ than in U_⊥, indicating that these white dwarfs likely do not originate from the local thin disc.

The AVRs of Holmberg et al. (2009) indicate that local thin disc stars are more dispersed in U than V, and we show in Fig. 3 that this expectation holds for the tangential velocity components U_⊥ and V_⊥ as well. Since the Holmberg et al. (2009) AVRs describe the full three-dimensional space velocities, we converted those relations to the expected AVRs of the projected velocities for the sky coordinate distribution of our 1.15–1.25 M_⊙ sample with 0.5–2.5 Gyr photometric ages. This enables us to properly compare the results of Holmberg et al. (2009) for main-sequence stars from the local thin disc to the dispersion values of the tangential velocity components inferred from Fig. 2 for our ultramassive white dwarf sample. The expected AVRs of the projected velocities are given by the square root of the diagonal elements of the expectation of the matrix |$\mathbf{ \boldsymbol{ \mathsf{ P}}} \mathbf{ \boldsymbol{ \mathsf{\sigma}}} ^2 \mathbf{ \boldsymbol{ \mathsf{ P}}}^{\mathrm T}$|⁠. Fig. 3 shows the AVRs from Holmberg et al. (2009) for the full three-dimensional space velocities for each of the U, V, and W components (as solid blue, orange, and green curves, respectively), compared to the computed AVRs expected from the projected velocities for our sample’s empirical sky distribution (dashed lines). We also analytically evaluated the expected AVRs of the projected velocities for an isotropic distribution (shown as the dotted curves in Fig. 3) and found them to be nearly the same as the AVRs determined using the empirical sky distribution of our sample.

Figure 3.

AVRs of projected and un-projected velocities. AVRs are shown for U, V, and W in blue, orange, and green, respectively. The solid curves show the AVRs for the full three-dimensional space velocities determined by Holmberg et al. (2009) for main-sequence stars. The dotted curves show the expected AVRs for the projected velocities assuming that the un-projected velocities follow the AVRs of Holmberg et al. (2009) and that the sky distribution is isotropic. The dashed curves also show the expected AVRs for the Holmberg et al. (2009) projected velocities but assuming that the sky distribution is that observed for 1.15–1.25 M_⊙ white dwarfs with 0.5–2.5 Gyr cooling ages.

Open in new tabDownload slide

From Fig. 3 we see that the expected dispersions for the projected velocities are typically notably smaller than the corresponding dispersion of the un-projected velocities, but that the relative size of the dispersion for V_⊥ compared to U_⊥ is similar to the relative sizes of σ_V and σ_U. Of particular note is that, throughout the age of the thin disc, the dispersion of V_⊥ is smaller than the dispersion of U_⊥. For the population of anomalously fast white dwarfs (mass bin 1.15–1.25 M_⊙ and photometric age bin 0.5–1.5 Gyr) with negative V_⊥, however, the dispersion of V_⊥ is actually larger than the dispersion of U_⊥. A population in equilibrium in the disc should have a dispersion ratio for V relative to U of σ_V/σ_U = 2/3 (Binney & Tremaine 2008), and the AVRs measured by Holmberg et al. (2009) for thin disc main-sequence stars have a dispersion ratio nearly equal to this equilibrium ratio (though actually slightly smaller). The dispersion of V_⊥ being larger than the dispersion of U_⊥ in the anomalous population thus indicates that this population is not in equilibrium in the disc and thus does not originate from the local thin (or thick) disc.

4 DISCUSSION

The presence of a population of anomalously fast-moving young, ultramassive white dwarfs on the Q branch was first discovered in Gaia DR2 observations by Cheng et al. (2019), who argued this population must experience an extra cooling delay on the Q branch of ∼6−8 Gyr that is not accounted for in current white dwarf cooling models. This proposed extra cooling delay scenario explained the large transverse velocities observed for the fast-moving population as a consequence of these white dwarfs being much older than their photometric ages indicate. The fast-movers in this scenario were assumed to originate from the local disc, where observed AVRs indicate that objects become more dispersed over time. As a check of the assumption of a local disc origin of the fast-mover population in the extra cooling delay scenario, Cheng et al. (2019) ran a test in which they modelled the velocity distribution of the fast-movers on the Q branch as a Gaussian and determined kinematic relations expected for a disc in equilibrium, notably finding a dispersion ratio σ_V/σ_U of approximately 2/3. Like Cheng et al. (2019), we also find a population of anomalously fast-moving young, ultramassive white dwarfs coincident with the Q branch. However, the empirical distributions of the transverse velocity components in Fig. 2 for Gaia EDR3 reveal kinematic features that are inconsistent with a local disc origin, in particular the distribution for V_⊥ and the dispersion ratio for V_⊥ and U_⊥.

Fig. 2 suggests that the distribution of V_⊥ for the anomalously fast ultramassive white dwarfs is actually too broad to be explained entirely by a cooling delay in (thin) disc white dwarfs. If white dwarfs with a true age of 10 Gyr produced through single stellar evolution in the Galactic thin disc follow the AVRs determined by Holmberg et al. (2009), then the expected dispersion for the projected velocity component V_⊥ in the tangent plane averaged over the sky distribution of the 1.15–1.25 M_⊙ white dwarfs with photometric ages 0.5–2.5 Gyr would be approximately 26 km s⁻¹ (see Fig. 3). From Fig. 2, we see that the anomalous 1.15–1.25 M_⊙ white dwarfs with 0.5–1.5 Gyr photometric ages have a dispersion for V_⊥ of ∼40 km/s when V_⊥ is negative (i.e. among white dwarfs that are moving in the opposite direction of Galactic rotation). The Galactic thin disc is not old enough to produce white dwarfs with such a high dispersion given the observed AVRs (Holmberg et al. 2009). While stars in the older thick disc have slightly higher dispersion (e.g. Sharma et al. 2014), the expected dispersion of ∼33 km s⁻¹ for V_⊥ is still too small to explain this dispersion through a cooling delay in a thick disc population.

Furthermore, as noted in Section 3.2 and seen in Fig. 2, we find a population of highly dispersed ultramassive white dwarfs coincident with the Q branch for which the dispersion in V_⊥ is larger than the dispersion in U_⊥. This is inconsistent with the dispersion ratio of a disc in equilibrium, for which the dispersion in V_⊥ should be smaller than the dispersion in U_⊥ with a dispersion ratio similar to that of the un-projected velocity components (see Fig. 3). Thus, the empirical dispersion ratio for V_⊥ and U_⊥ suggests that these white dwarfs do not originate from the Milky Way disc.

The true origin of this population, however, is unclear. These white dwarfs may have originated from elsewhere in the Galaxy such as from the halo. The halo is older than the Galactic disc (e.g. Gallart et al. 2019) and halo stars are typically more dispersed (Gaia Collaboration 2018b). This scenario still requires a mechanism by which photometrically young white dwarfs acquire velocities faster than their photometric age would indicate. This mechanism could be an extra cooling delay as proposed by Cheng et al. (2019) of nearly 10 Gyr, so that only white dwarfs with photometric ages of 0.5–1.5 Gyr are affected. Alternatively, these objects could originate from a dynamically distinct, hot yet young population.

Another possible explanation of the anomalous fast-movers is that they are the merger remnant of the inner binary of a hierarchical triple system whose outer companion was ejected. Observations indicate that triple systems are common, constituting about 10 per cent of multiples with F and G dwarf primaries (Eggleton & Tokovinin 2008; Raghavan et al. 2010; Tokovinin 2014; Moe & Di Stefano 2017) and much larger fractions for higher mass B- and O-type primaries (Sana et al. 2014; Moe & Di Stefano 2017), and population synthesis simulations indicate that stellar interactions occur in the majority of triple systems (Toonen et al. 2020; Stegmann, Antonini & Moe 2022). The evolution of triple systems is significantly more complicated than that of isolated binaries due to the combination of stellar evolution, stellar interactions, and three-body dynamics (see e.g. Toonen, Hamers & Portegies Zwart 2016 for a review). Triple systems tend to be hierarchical, with an inner binary and a distant companion (Hut & Bahcall 1983). Such systems have a large variety of possible evolution channels (Toonen et al. 2020; Stegmann et al. 2022; Toonen, Boekholt & Portegies Zwart 2022) and can result in mergers and ejections (Mardling & Aarseth 2001; Perets & Kratter 2012; Glanz & Perets 2021; Toonen et al. 2022; Grishin & Perets 2022; Hamers, Glanz & Neunteufel 2022a; Hamers et al. 2022b). The secular gravitational effects between the outer companion and inner binary can drive the inner binary to high eccentricity in Zeipel–Lidov–Kozai oscillations (von Zeipel 1910; Kozai 1962; Lidov 1962; Naoz 2016; Ito & Ohtsuka 2019), which can enhance mergers of compact objects (e.g. Blaes, Lee & Socrates 2002; Thompson 2011; Hamers et al. 2013; Antognini et al. 2014; Antonini, Toonen & Hamers 2017; Liu & Lai 2017, 2018; Silsbee & Tremaine 2017; Hamers et al. 2018; Hoang et al. 2018; Randall & Xianyu 2018a, b; Toonen, Perets & Hamers 2018; Fragione & Loeb 2019; Stegmann et al. 2022). For example, dynamical effects of the outer companion can enhance the merger rate of inner compact double white dwarf binaries relative to isolated binary systems (Thompson 2011). Destabilized stellar triple systems in the field can eject a component with runaway speeds of tens of km s⁻¹ (Toonen et al. 2022).

Unstable triple systems in star clusters can also result in the ejection of the tertiary perturber and escape from the cluster (Mardling & Aarseth 2001). From their N-body simulations of open clusters, Mardling & Aarseth (2001) describe a particular event in which an initially stable triple system became unstable and the tertiary companion was ejected; the resulting binary and single star both escaped the cluster with respective terminal velocities of 20.6 and 65.7 km s⁻¹. Star clusters are known to harbour exotica such as blue stragglers, which can be formed through the dynamical evolution of stellar triple systems (Iben & Tutukov 1999) and have been found in both open clusters (Ahumada & Lapasset 2005; de Marchi et al. 2006) and globular clusters (Ferraro et al. 2012; Parada et al. 2016). Furthermore, cluster members can display different kinematic relations than what is expected for a disc in equilibrium, and many nearby clusters are moving at speeds of tens of km s⁻¹ relative to the local standard of rest (Gaia Collaboration 2018a; Soubiran et al. 2018, 2019; Lodieu et al. 2019). The kinematically anomalous population we observe could thus conceivably originate from old open clusters or globular clusters that are close to the Galactic disc. While open clusters are typically younger than the members of the anomalous population, some open clusters are old enough to be a source of these white dwarfs, and those old open clusters are also typically rich (Gaia Collaboration 2018a; Soubiran et al. 2018, 2019). As open clusters are only loosely gravitationally bound, stars can also readily escape from them; for example, observations indicate that the Pleiades has lost ∼20 per cent of its mass over the past 100 Myr (Heyl, Caiazzo & Richer 2022). Though globular clusters are more tightly bound than open clusters and thus less prone to the escape of cluster members, they are also typically much more populated (Gaia Collaboration 2018a) and only a small fraction of escapees is needed to explain the anomalous population that we observe in this work.

In order to explain the clustering in photometric age in this scenario, the fast-moving population must be produced through a singular event rather than some continuous process. This could be cluster evaporation in the case of an open cluster origin, or it could be some transient event that disrupted a globular cluster. This transient event may relate to the burst of star formation in the local Galactic disc that peaked around 2–3 Gyr ago and continued until ∼1 Gyr ago (Mor et al. 2019). A travel time of around 1 Gyr from the location of origin after the disruption of the cluster would then produce the observed clustering in age of the fast-moving population.

5 CONCLUSIONS

As a follow-up to our previous work on the distribution of photometric cooling ages of massive white dwarfs in Gaia EDR3 (Paper I), we analysed the kinematics of the transverse motion of 0.95–1.25 M_⊙ Gaia EDR3 white dwarfs. We find that the population of anomalously fast ultramassive white dwarfs on the Q branch reported by Cheng et al. (2019) for Gaia DR2 is still present in Gaia EDR3. These white dwarfs appear photometrically young according to single stellar evolution models, but are moving faster than expected for stars of that age according to AVRs observed for the local thin disc. Using white dwarf cooling models for which the photometric age distributions are consistent with the expectation from the star formation history (Paper I) to infer the masses and ages of our sample of white dwarfs and sorting our full sample into mass and age bins, we find that this population of fast-movers is concentrated to masses of 1.15–1.25 M_⊙ and photometric cooling ages of 0.5–1.5 Gyr.

Our analysis of the distributions of the individual components of the transverse velocity reveals that among the white dwarfs in the mass range 1.15–1.25 M_⊙ and age range 0.5–1.5 Gyr, there is a population of white dwarfs lagging the local Galactic rotation that is too dispersed in V_⊥ to be explained solely by a cooling delay in white dwarfs born in the local thin disc. According to the AVRs determined from observations of main-sequence stars in the local thin disc (Holmberg et al. 2009), it would take longer than the age of the thin disc (10 Gyr) for disc heating to produce a dispersion large enough to explain this anomalous population of white dwarfs. Furthermore, we find this population to be more dispersed in V_⊥ than in U_⊥, which suggests that this population does not originate from the local disc. Some potential explanations of this population include a halo origin, in conjunction with an extra cooling delay or from a dynamically distinct population, or that they are produced through triple system dynamics such as the merger of the inner binary of a hierarchical triple system whose tertiary companion was ejected resulting in a large velocity for the binary that ultimately merges to form the white dwarf. However, the precise origin of this population of young ultramassive white dwarfs is an open question for future work.

ACKNOWLEDGEMENTS

This work has been supported by the Natural Sciences and Engineering Research Council of Canada through the Discovery Grants program and Compute Canada. IC is a Sherman Fairchild Fellow at Caltech and thanks the Burke Institute at Caltech for supporting her research.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We used the ‘Synthetic Colors and Evolutionary Sequences of Hydrogen- and Helium-Atmosphere White Dwarfs’ website at http://www.astro.umontreal.ca/~bergeron/CoolingModels/.

DATA AVAILABILITY

The data used in this paper are available through TAP Vizier and the Gaia archive. We constructed the 200-pc white dwarf catalogue from https://warwick.ac.uk/fac/sci/physics/research/astro/research/catalogues/gaiaedr3_wd_main.fits.gz. The corresponding author can also provide software to perform the analysis presented in this paper.

Footnotes

REFERENCES

Author notes