A Cooling Anomaly of High-mass White Dwarfs

Gentile Fusillo et al. (2019) have compiled a catalog of Gaia DR2 white dwarfs based on the G -band absolute magnitude, Gaia color index, and some quality cuts. To select white dwarfs with high-precision measurements, we further apply the following quality cuts:

$$\begin{eqnarray}&&{\sigma }_{{G}_{\mathrm{BP}}-{G}_{\mathrm{RP}}}\lt 0.10,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\sigma }_{\mu }/\varpi \lt 2\,\mathrm{km}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\mathtt{parallax}}{\mathtt{\_}}{\mathtt{over}}{\mathtt{\_}}{\mathtt{error}}\gt 8,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\mathtt{a}}{\mathtt{s}}{\mathtt{t}}{\mathtt{r}}{\mathtt{o}}{\mathtt{m}}{\mathtt{e}}{\mathtt{t}}{\mathtt{r}}{\mathtt{i}}{\mathtt{c}}{\rm{\_}}{\mathtt{e}}{\mathtt{x}}{\mathtt{c}}{\mathtt{e}}{\mathtt{s}}{\mathtt{s}}{\rm{\_}}{\mathtt{n}}{\mathtt{o}}{\mathtt{i}}{\mathtt{s}}{\mathtt{e}}\lt 1.5,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\mathtt{phot}}{\rm{\_}}{\mathtt{bp}}{\rm{\_}}{\mathtt{rp}}{\rm{\_}}{\mathtt{excess}}{\rm{\_}}{\mathtt{factor}}\lt 1.4,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&1/\varpi \lt 250\,{\rm{p}}{\rm{c}},\end{eqnarray} \tag{ 6 }$$

where the color error ${\sigma }_{{G}_{\mathrm{BP}}-{G}_{\mathrm{RP}}}$ is the combined photometric errors in G_BP and G_RP bands, the proper motion error σ_μ is the combined error originating from its two components, and ϖ is the parallax from Gaia DR2. These cuts are designed to balance data quality and sample size. They do not introduce explicit kinematic biases, which is necessary for our analysis below. While the main sample used in our study uses white dwarfs within 250 pc, we will occasionally use a subsample of white dwarfs located within 150 pc to clearly show number density enhancements.

Our analysis requires white dwarf absolute magnitude, color index, and the two components of transverse velocity ${{\boldsymbol{v}}}_{{\rm{T}}}$ = (${v}_{{\rm{L}}}$, ${v}_{{\rm{B}}}$) in Galactic longitude and latitude directions. Except for the color index G_BP–G_RP, which is directly read from the bp_rp column in Gaia DR2, we derive the other quantities in the following way:

$$\begin{eqnarray}&&{M}_{G}=G+5\mathrm{log}(\varpi /\mathrm{mas})-10,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{v}_{{\rm{L}}}=\displaystyle \frac{{\mu }_{{\rm{L}}}-(A\cos 2l+B)\cos b}{\varpi },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{v}_{{\rm{B}}}=\displaystyle \frac{{\mu }_{{\rm{B}}}+A\sin 2l\sin b\cos b}{\varpi },\end{eqnarray} \tag{ 9 }$$

where G and ϖ/mas are read from Gaia DR2 columns phot_g_mean_mag and parallax, μ_L and μ_B are converted from columns ra, dec, pmra, and pmdec with the coordinate conversion function in the astropy package (Astropy Collaboration et al. 2013, 2018), and A and B are the Oort constants taken from Bovy (2017). We do not correct for extinction because within the distance cut, extinction is in general tiny and there is no accurate estimate for it. To avoid the influence of hyper-velocity white dwarfs, we further impose a velocity cut:

$$\begin{eqnarray}&&{v}_{{\rm{T}}}=\sqrt{{v}_{{\rm{L}}}^{2}+{v}_{{\rm{B}}}^{2}}\lt 250\,\mathrm{km}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 10 }$$

We point out that Gaia does not provide any useful radial velocity information for white dwarfs as they have no spectral lines in the 845–872 nm wavelength range of Gaia's spectrometer (Gaia Collaboration et al. 2016).

We then derive white dwarf photometric isochrone ages and masses from the H–R diagram coordinates:

$$\begin{eqnarray}&&({G}_{\mathrm{BP}}\mbox{--}{G}_{\mathrm{RP}},{M}_{G})\to ({\tau }_{\mathrm{phot}},{m}_{\mathrm{WD}}),\end{eqnarray} \tag{ 11 }$$

based on a single-star evolution scenario and white dwarf cooling models.⁴ We estimate the main-sequence (MS) ages with an initial–final mass relation (Cummings et al. 2018) and the relation between pre-WD time and main-sequence mass from Choi et al. (2016) for non-rotating, solar-metallicity stars. For high-mass white dwarfs, the pre-WD ages are negligible. As for white dwarf cooling, we use a table of synthetic colors for pure-hydrogen atmosphere (Holberg & Bergeron 2006; Kowalski & Saumon 2006; Tremblay et al. 2011) and a grid of cooling tracks for C/O-core white dwarfs with "thick" hydrogen layers (Fontaine et al. 2001).⁵ In order to convert any H–R diagram coordinate into τ_phot and m_WD, we linearly interpolate between grid points. Stellar models show that in the single-star-evolution scenario, white dwarfs with a mass higher than about 1.05–1.10 M_⊙ have oxygen+neon (O/Ne) cores (e.g., Siess 2007; Lauffer et al. 2018). So, we combine the cooling tracks of m_WD ≤ 1.05 M_⊙ C/O white dwarfs with the four cooling tracks of m_WD ≥ 1.10 M_⊙ O/Ne white dwarfs (Camisassa et al. 2019).

The O/Ne white dwarf model only gives slightly lower mass estimates than the C/O white dwarf model (e.g., 1.08–1.23 M_⊙ in the combined model corresponds to 1.10–1.28 M_⊙ in the C/O-only model), and their estimates of the photometric ages τ_phot are similar for the white dwarfs we are interested in (τ_phot < 3.5 Gyr). Switching between thick-hydrogen, thin-hydrogen, and helium atmosphere (Bergeron et al. 2011) models does not significantly change the photometric-age estimate of our sample either.

In Figure 1 we show the selected white dwarfs on the H–R diagram. In the top right panel, we use the 150 pc sample to show the density distribution with a higher contrast. The Q branch is a factor-two enhancement at around $-0.4\,\lt {G}_{\mathrm{BP}}\mbox{--}{G}_{\mathrm{RP}}\lt 0.2$ and M_G = 13. In the main panel, we show our main sample within 250 pc, color-coded by their transverse velocities v_T with respect to the local standard of rest. We adopt (U_⊙, V_⊙, W_⊙) = (11, 12, 7) km s⁻¹ (Schönrich et al. 2010) to correct for the solar reflex motion. We emphasize the fast white dwarfs (v_T > 70 km s⁻¹) in Figure 1 with larger dots: they are very likely thick-disk stars. We also plot a grid of photometric age τ_phot and mass m_WD derived from the combined O/Ne-core and C/O-core white dwarf cooling model. Cooling tracks are the curves with constant m_WD. White dwarfs with different birth times form a "white dwarf cooling flow" on the H–R diagram as they move along their cooling tracks.

Zoom In Zoom Out Reset image size

We focus on the mass range where the Q branch is most prominent. To maximize sample size and minimize the contamination from standard-mass helium-atmosphere white dwarfs (the B branch), we impose the following photometric age and mass cuts:

$$\begin{eqnarray}&&0.1\,\mathrm{Gyr}\lt \,{\tau }_{\mathrm{phot}}\lt 3.5\,\mathrm{Gyr},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&1.08\,{M}_{\odot }\lt \,{m}_{\mathrm{WD}}\lt 1.23\,{M}_{\odot },\end{eqnarray} \tag{ 13 }$$

where m_WD is derived from the combined cooling model for O/Ne and C/O white dwarfs. This mass range corresponds to 1.10–1.28 M_⊙ in the C/O-only cooling model. In total, 1070 white dwarfs are selected by these criteria.⁶ In this region, the Q branch divides the white dwarf cooling flow into three segments: the early, Q-branch, and late segments, as shown in Figure 1. We define the Q-branch segment by

$$\begin{eqnarray}&&| {M}_{G}-1.2\times {\mathtt{bp}}\_{\mathtt{rp}}-13.2| \lt 0.2\end{eqnarray} \tag{ 14 }$$

in addition to the previous photometric-age and mass cuts.

As argued by Tremblay et al. (2019), an enhancement not aligned with mass or age grid, such as the Q branch, should be produced by a slowing down (and therefore a delay) of white dwarf cooling. Such a cooling delay creates a "traffic jam" in the white dwarf flow, and the Q branch is a snapshot of this traffic jam. Is crystallization alone enough to explain the cooling delay on the Q branch? If it is, then the distribution of photometric age τ_phot derived from a cooling model including crystallization effects should no longer carry signatures of the Q branch. However, observations lead to the antithesis. In Figure 2 we show the distribution of τ_phot in three mass ranges: there is a mass-dependent enhancement tracing the Q branch, which is consistent with the observation by Tremblay et al. (2019) that the pile-up is higher and narrower than what the standard cooling model predicts. Evolutionary delays from binary interactions or a peak in star formation rate cannot explain this mass-dependent enhancement either. Therefore, an extra cooling delay in addition to crystallization effects (latent heat and phase separation) must exist.

Zoom In Zoom Out Reset image size

Observations show that the velocity dispersion of disk stars in the Milky Way is related to the stellar age τ: older stars have higher velocity dispersion than younger stars (e.g., Holmberg et al. 2009). So, the transverse velocity ${{\boldsymbol{v}}}_{{\rm{T}}}$ derived from Gaia DR2 can be used as a "dynamical" indicator of the true stellar age τ. For the Milky Way thin disk, the dispersion of transverse velocity approximately follows a power law increasing from about 25 km s⁻¹ at 1.5 Gyr to 55 km s⁻¹ at around 6–8 Gyr (e.g., Holmberg et al. 2009); for the thick disk, the dispersion is about 65 km s⁻¹ (e.g., Sharma et al. 2014). Given this age–velocity-dispersion relation (AVR), we observe two anomalous things in the velocity distribution of the Q branch:

• 1.  
  There is a strong excess of white dwarfs with v_T > 70 km s⁻¹ in the Q-branch segment, as shown in Figure 1. According to the age–velocity-dispersion relation (AVR) mentioned above, these fast white dwarfs should be old stars. Given that the photometric age on the Q branch is only 0.5–2 Gyr, these white dwarfs must have experienced an extra cooling delay for several billion years. In the left panel of Figure 3 we show that the excess of fast white dwarfs in the Q-branch segment is clear for m_WD > 1.05 M_⊙; in the right panel we show that this excess is observed for a variety of velocity cuts.
• 2.  
  The fraction of fast white dwarfs in the late segment is lower than that in the Q-branch segment. This is anomalous, because white dwarfs in the late segment should be older than those in the Q-branch segment, as long as all white dwarfs follow the same cooling law. The only way to create such a reverse of fraction is to have more than one white dwarf population with distinct cooling behaviors.

Zoom In Zoom Out Reset image size

The simplest scenario that can explain both the number enhancement and velocity anomaly on the Q branch is to have an extra-delayed population of white dwarfs in addition to a normal-cooling population. This scenario requires only two free parameters:

• 1.  
  The fraction f_extra of the extra-delayed population, and
• 2.  
  The length t_extra of the extra cooling delay on the branch.

In Figure 4 we illustrate this scenario by showing the normal-cooling and extra-delayed populations on the H–R diagram. Before the Q branch, the extra-delayed population has no difference from the normal-cooling population. On the branch, the extra-delayed population has a slower cooling rate, which causes two effects: (1) its members pile up there, creating a number-density enhancement, and (2) the photometric ages τ_phot of its members start to seem younger than their true ages τ, creating an age discrepancy. After the branch, the number-density enhancement disappears, but the age discrepancy remains. A detailed parameterization of this scenario is presented in Appendix A. To create the observed reversal of fast white dwarf fraction in the Q-branch and late segments, the extra-delayed population must have a high number-density contrast between these two regions, which requires that the population fraction f_extra be small and the delay time t_extra long.

Simulations of binary evolution show that double-WD merger products may account for a considerable fraction of high-mass white dwarfs (e.g., Toonen et al. 2017; Temmink et al. 2019). These merger products also have a discrepancy between their true ages and photometric ages due to binary evolution. Therefore, in order to use the velocity distribution to quantitatively constrain the extra cooling delay, the merger population must be modeled simultaneously. Constraining the merger fraction is also of great interest as its value is still a matter of debate (e.g., Giammichele et al. 2012; Wegg & Phinney 2012; Rebassa-Mansergas et al. 2015; Tremblay et al. 2016; Maoz et al. 2018). Therefore, we include the double-WD merger products in our model and set their fraction as a free parameter.

In our model, we consider two evolutionary delays: the extra cooling delay, and the merger delay. Accordingly, we consider three populations of white dwarfs with different combinations of the two delays:

• 1.  
  A generic population of singly evolved white dwarfs that follows normal cooling, denoted by "s;"
• 2.  
  A double-WD merger population⁷ with systematic age offsets due to the merger delay and with a normal cooling, denoted by "m;"
• 3.  
  A population with the extra cooling delay, denoted by "extra."

Their delay scenarios are listed in Table 1. For simplicity, we only explore the two extreme situations for the extra-delayed population, where

• 1.  
  Setup 1: all members of the extra-delayed population also have the merger delay;
• 2.  
  Setup 2: no members of the extra-delayed population have the merger delay.

The distribution function p(${\boldsymbol{y}}$) of observables ${\boldsymbol{y}}$ for all white dwarfs can be written as a weighted average of the distribution for each population:

$$\begin{eqnarray}&&p({\boldsymbol{y}})={f}_{{\rm{s}}}\,{p}_{{\rm{s}}}({\boldsymbol{y}})+{f}_{{\rm{m}}}\,{p}_{{\rm{m}}}({\boldsymbol{y}})+{f}_{\mathrm{extra}}\,{p}_{\mathrm{extra}}({\boldsymbol{y}}),\end{eqnarray} \tag{ 15 }$$

where the weight f denotes the fraction of each population, satisfying f_s + f_m + f_extra = 1 .

Our goal is to use observations to constrain two independent population fractions and the delay time of the extra cooling delay:

$$\begin{eqnarray}&&{f}_{{\rm{m}}},{f}_{\mathrm{extra}},\ \mathrm{and}\ {t}_{\mathrm{extra}},\end{eqnarray} \tag{ 16 }$$

the last of which is encoded in the distribution p_extra(y). We have two sets of observables y: the transverse velocities v_T, and the photometric ages τ_phot. They are connected by the AVR $p({\boldsymbol{v}}| \tau )$ and the delay scenario of each population (listed in Table 1). The delay

$$\begin{eqnarray}&&{\rm{\Delta }}t\equiv \tau -{\tau }_{\mathrm{phot}}\,\end{eqnarray} \tag{ 17 }$$

includes contributions from the extra-cooling and/or the merger delays. We build a Bayesian model based on Equation (15) to constrain the aforementioned parameters. Our model is similar to that of Wegg & Phinney (2012), but we include the extra-delayed population and use a much larger sample. In addition, to avoid the need for modeling selection effects, we derive our constraints from the velocity distribution conditioned on observables other than velocity:

$$\begin{eqnarray}&&p({\boldsymbol{v}}| \mathrm{other}\,\mathrm{observables})=p({{\boldsymbol{v}}}_{{\rm{T}}}| {\tau }_{\mathrm{phot}},{m}_{\mathrm{WD}},l,b).\end{eqnarray} \tag{ 18 }$$

The details of this statistical technique and the Bayesian framework of our model are shown in Appendix B.

The free parameters in our model include the population fractions f_m and f_extra, the extra delay time t_extra, parameters for the AVR, and solar motion. Although constraints on the AVR and solar motion already exist, treating them also as free parameters can avoid potential systematic errors, and the comparison of our best-fitting values with the existing values allows us to check the validity of our method.

Below, we list the main assumptions and simplifications in our model:

• 1.  
  We assume that upon entering the Q-branch segment, all members of the extra-delayed population suddenly slow down their cooling by a constant factor, and upon leaving the branch, the cooling rates suddenly resume, so that this extra cooling delay can be parameterized by just its length t_extra and population fraction f_extra (see Appendix A). The resulting delay-time distribution is described in Section 4.3.
• 2.  
  The velocity distribution of white dwarfs is a superposition of 3D Gaussian distributions as a function of age τ, i.e., $p({\boldsymbol{v}}| \tau )={ \mathcal N }({{\boldsymbol{v}}}_{0}(\tau ),{\boldsymbol{\Sigma }}(\tau ))$. The details of this Gaussian velocity model are shown in Appendix C.
• 3.  
  The true-age distribution of high-mass white dwarfs within 250 pc is uniform up to 11 Gyr, i.e., τ ∼ U [0, 11 Gyr].
• 4.  
  For the double-WD merger products, we follow Wegg & Phinney (2012) and assume that the resulting white dwarf is reheated enough that its cooling age after the merger is almost equal to the photometric cooling age. We also assume a fixed delay-time distribution for double-WD mergers (see Section 4.3) and a parameterization of the AVR (see Section 4.4).

The three white dwarf populations in our model are defined by their different delay signatures Δt, which concern the extra cooling delay Δt_extra and double-WD merger delay Δt_merger. The delay scenario of each population in each segment is listed in Table 1.

The extra cooling delay Δt_extra is built up on the Q branch. We adopt a uniform distribution Δt_extra ∼ U [0, t_extra] of this delay for white dwarfs in the Q-branch segment and a constant value in the late segment. Note that Δt_extra is a random variable with a probability distribution, whereas t_extra, as a model parameter to be constrained, is the upper limit of Δt_extra. In the Q-branch segment, we do not further distinguish if a white dwarf has just started or is about to complete their extra cooling delay, because the uncertainty of H–R diagram coordinate due to different atmosphere types and astrometric/photometric error is comparable to the width of the Q branch on the H–R diagram. In this case, a uniform distribution is a good and efficient approximation for Δt_extra.

The double-WD merger delay Δt_merger originates from the binary evolution before the merger. We refer to binary population synthesis results (e.g., Toonen et al. 2014) and approximate the delay by

$$\begin{eqnarray}&&p({\rm{\Delta }}{t}_{\mathrm{merger}})\propto {\rm{\Delta }}{t}_{\mathrm{merger}}^{-0.7}\end{eqnarray} \tag{ 19 }$$

for Δt_merger > 0.5 Gyr and zero for smaller Δt_merger. Unlike the extra cooling delay, we do not set any free parameter for this merger-delay distribution.

We define the U, V, W axes as pointing toward (l = 0°, b = 0°), (l = 90°, b = 0°), and (b = 90°), respectively, and assume that the main axes of the Gaussian velocity distribution are aligned with these directions with dispersion σ_U(τ), σ_V(τ), and σ_W(τ). Observations show that the AVR in each direction can be fit by a shifted power law. The power index of the in-disk components are around 0.35 and that of the W component is around 0.5 (e.g., Holmberg et al. 2009; Sharma et al. 2014). For old stars including thick-disk members, the AVR is still a matter of debate (e.g., Yu & Liu 2018; Mackereth et al. 2019). So, in each direction, we use a shifted power law to parameterize the AVR of the younger, thin-disk stars:

$$\begin{eqnarray}&&\sigma (\tau )={\sigma }^{\tau =0}+{\rm{\Delta }}{\sigma }^{0\to 4}\times {\left(\displaystyle \frac{\tau }{4}\right)}^{\beta },\ \ \ \ \tau \lt 7\,\mathrm{Gyr},\end{eqnarray} \tag{ 20 }$$

and we use a constant value σ^thick to represent the velocity dispersion of stars older than 10 Gyr (thick-disk stars); in between 7 and 10 Gyr, we linearly interpolate the values from the two ends to reflect the increasing fraction of thick-disk stars. The shape of the AVR with our parameterization is shown in Figure 5. The ratio of the two in-disk components σ_V and σ_U should be a constant for a local sample (e.g., Binney & Tremaine 2008), so we set σ_V(τ) = k σ_U(τ). As the assumption for the velocity distribution to be Gaussian gradually breaks down when σ_U increases, we allow the ratio k to be different for the thin and thick disks. Thus, we use in total 10 parameters to model the anisotropic AVR: two initial velocity dispersions ${\sigma }_{{\text{}}U,W}^{\tau =0}$, two dispersion increases ${\rm{\Delta }}{\sigma }_{{\text{}}U,W}^{0\to 4}$ between 0 and 4 Gyr, two power indices β_{U, W}, two thick-disk dispersions ${\sigma }_{{\text{}}U,W}^{\mathrm{thick}}$, and two in-disk dispersion ratios k^thin and k^thick. The best-fitting values of these parameters can be checked against existing estimates presented in the literature.

Zoom In Zoom Out Reset image size

In Figure 6 we present the constraints we obtain for the parameters of interest: f_extra, t_extra, and f_m. We find that the extra-delayed population fraction is

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{extra}} & = & {6.4}_{-1.5}^{+2.3} \% \,\,(\mathrm{setup}\ 1)\\ & = & {9.2}_{-2.7}^{+4.4} \% \,\,(\mathrm{setup}\ 2),\end{array}\end{eqnarray} \tag{ 21 }$$

and the length of the extra cooling delay is

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{extra}} & \geqslant & 8\,\mathrm{Gyr}\ \,(\mathrm{setup}\ 1)\\ & \geqslant & 10\,\mathrm{Gyr}\,\,(\mathrm{setup}\ 2).\end{array}\end{eqnarray} \tag{ 22 }$$

These constrains confirm our qualitative conclusion that f_extra is small and t_extra is long in Section 3.3. We point out that the difference of t_extra in the two setups is exactly where the peak of the merger-delay distribution is located (2 Gyr), which is expected. The lower limit for t_extra slightly depends on the parameterization of the AVR: if we adopt a younger thick-disk age (7–11 Gyr instead of 9–11 Gyr, Mackereth et al. 2019), this lower limit will also decrease by about 1–2 Gyr. ${t}_{{\rm{e}}{\rm{x}}{\rm{t}}{\rm{r}}{\rm{a}}}$ may be overestimated for the fact that at a high level of dispersion, the velocity distribution is not exactly Gaussian. But it is unlikely that we overestimate ${t}_{{\rm{e}}{\rm{x}}{\rm{t}}{\rm{r}}{\rm{a}}}$ too much, because we set the thick-disk dispersion as a free parameter. We have also checked that reasonable variations of the input delay-time distribution of the mergers (e.g., Toonen et al. 2012) do not change these two constraints significantly.

Zoom In Zoom Out Reset image size

The fraction of merger products without the extra cooling delay is found to be ${f}_{{\rm{m}}}={15}_{-5}^{+6} \% \ \mathrm{and}\ {20}_{-5}^{+6} \% \,$ (setups 1 and 2). Therefore, the total fraction of double-WD merger products is

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{extra}}+{f}_{{\rm{m}}} & = & {22}_{-5}^{+7} \% \,(\mathrm{setup}\ 1)\\ {f}_{{\rm{m}}} & = & {20}_{-5}^{+6} \% \,(\mathrm{setup}\ 2)\end{array}\end{eqnarray} \tag{ 23 }$$

among 1.08–1.23 M_⊙ white dwarfs. This total fraction is mainly constrained by the fast white dwarfs in the early segment (where the two setups do not differ from each other), so the constraints on this fraction under setups 1 and 2 are similar. A more detailed analysis of the merger products among high-mass white dwarfs is presented in Cheng et al. (2019).

Finally, we calculate the contribution of the extra-delayed population in the Q-branch segment according to the above fractions:

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{extra}} & = & \displaystyle \frac{{f}_{\mathrm{extra}}({t}_{\mathrm{extra}}/{\rm{\Delta }}{t}_{\mathrm{branch}}+1)}{1+{f}_{\mathrm{extra}}{t}_{\mathrm{extra}}/{\rm{\Delta }}{t}_{\mathrm{branch}}}\\ & = & (47\pm 8) \% ,\end{array}\end{eqnarray} \tag{ 24 }$$

where Δt_branch is the width of the Q branch (see Appendix A). This fraction, derived purely from the velocity information, is consistent with the factor ∼2 estimate obtained directly from the number-density enhancement in Figure 2, which confirms that our extra-cooling-delay scenario is a good phenomenological model for the Q branch.

To further validate our model and results, we first check our constraints on the nuisance parameters. For the solar motion we obtain

$$\begin{eqnarray}&&{\left(U,V,W\right)}_{\odot }=(10.3\pm 1.0,7.3\pm 1.0,6.7\pm 0.5)\,\mathrm{km}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 25 }$$

which is consistent with the measurement of Rowell & Kilic (2019) based on mainly standard-mass white dwarfs. Our values of U_⊙ and W_⊙ are also consistent with the results of Schönrich et al. (2010). The discrepancy of the V_⊙ measurement comes from the different treatments of the asymmetric drift and is beyond the scope of this paper.

Figure 5 shows our constraint on the white-dwarf AVR, which is consistent with the AVR of thin- and thick-disk MS stars (Holmberg et al. 2009; Sharma et al. 2014). Removing either the extra-delayed population or the merger population leads to unreasonably higher AVRs. Before Gaia DR2 came out, Anguiano et al. (2017) reported an unexpectedly high AVR for young white dwarfs (their Figures 21 and 22), without considering the extra cooling delay or merger delay in their age estimate. This unreasonably high AVR is exactly what we see when we remove the extra-delayed and/or merger population from our model. Thus, we verify that both the extra cooling delay and the merger delay are necessary.

To check the goodness of fit, we compare the observed and modeled velocity distributions in Figure 7. Our best-fitting models (in both setups 1 and 2) provide good characterizations of the observed velocity distribution in all the early, Q-branch, and late segments. Adopting a different star formation history introduces no significant changes to our results. We test both a linearly decreasing star formation rate with a five-time higher star formation rate in the past, and a star formation history with a bump at 2.5 Gyr ago (e.g., Mor et al. 2019), and find that the changes in best-fitting values are smaller than their uncertainties. This insensitivity to the assumed star formation history is expected because our model mainly uses the velocity information.

Zoom In Zoom Out Reset image size

To further argue for our extra-delayed scenario against other explanations of the velocity anomaly, such as an accretion event of the Milky Way, we run a simple test where the velocities of fast white dwarfs on the Q branch are parameterized by only one Gaussian distribution. We find that the mean of the U and W components are consistent with zero, and the mean of V is −50 ± 6 km s⁻¹. Moreover, the U component has a dispersion of 60 ± 6 km s⁻¹, and the ratio between the U and V dispersion is 0.60 ± 0.08. All of these values satisfy the relations for a disk in equilibrium (e.g., Binney & Tremaine 2008): asymmetric drift $\bar{V}=-{\sigma }_{{\text{}}U}^{2}/(80\,\mathrm{km}\,{{\rm{s}}}^{-1})$ and dispersion ratio σ_V/σ_U = 0.67. It is unlikely for accreted stars to exactly reproduce the disk kinematics.

The Q branch is named after the presence of DQ-type white dwarfs (Gaia Collaboration et al. 2018a). To explore this dimension, we cross-match our white dwarf sample with the Montreal white dwarf database, MWDD (Dufour et al. 2017).⁸ We note that most high-mass DQs are concentrated on the branch (Figure 9) and the fraction of fast DQs on the branch is very high (Table 2). Therefore, all of these Q-branch DQs are likely to belong to the extra-delayed population. However, not all extra-delayed white dwarfs are DQs. We estimate the fraction of DQs in the extra-delayed population to be

$$\begin{eqnarray}&&{F}_{\mathrm{DQ}}=19/(76{F}_{\mathrm{extra}})=(53\pm 16) \% ,\end{eqnarray} \tag{ 26 }$$

based on the total number of DQs and DAs in Table 2. Changing the distance limit of the sample does not influence this result much. It remains unclear but is of further interest to investigate the reason why half of the extra-delayed white dwarfs are DQs while the other half are DAs.

Zoom In Zoom Out Reset image size

The DQs on the Q-branch are abnormal because the convection zone in a normal white dwarf with similar temperature is not deep enough to dredge up carbon (Dufour et al. 2005). In a similar way, the hot-DQ white dwarfs discovered by Dufour et al. (2007) are also abnormal. In Figure 9 we show the distributions of the Q-branch DQs and hot DQs on the H–R diagram. Below, we argue that although these two groups of DQs are observationally different, they may be related through an evolutionary relation.

The hot DQs and Q-branch DQs appear to be different in some aspects. Hot-DQ white dwarfs are characterized by the high temperature (>18,000 K), highly carbon-dominant atmosphere (Williams et al. 2013), high rate of having a magnetic field (Dufour et al. 2010, 2013), being variable (e.g., Dufour et al. 2009, 2011; Dunlap et al. 2010; Williams et al. 2016), and rarity (e.g., Dufour et al. 2008). In contrast, the Q-branch DQ white dwarfs are concentrated on the Q branch, have helium-dominant atmospheres with 10⁻⁴–10⁻¹ carbon (Kepler et al. 2015, 2016; Coutu et al. 2019), and have undetectable or no magnetic field (see Figure 9; a caveat for the magnetic field is that most hot DQs have been examined with high-resolution spectroscopy, so their magnetic fields are more likely to be found). As for kinematics, hot DQs are mildly faster than normal white dwarfs, which is an indication of being merger products (Dunlap & Clemens 2015), whereas Q-branch DQs are much faster, which needs the long extra cooling delay to explain. Dunlap & Clemens (2015) discussed one strange hot DQ (SDSS J115305.47+005645.8 or J1153) with a very high proper motion. We note that J1153 has not been reported to have magnetic field or variability and lies on the Q branch (Figure 9), which means that J1153 can be classified as a Q-branch DQ.

We now turn to the similarities between Q-branch DQs and hot DQs: both of them have high masses, and both of them might have a merger origin. These two similarities raise a serious question: are the Q-branch DQs evolved from the hot DQs? We use number counts to explore this possibility. Hot DQs are rare; based on a spectroscopically verified white dwarfs sample (as shown in Figure 9), we find that the fraction of hot DQs is 8/203 = 4.0 ± 1.4% in the region earlier than the Q branch (τ_phot < 0.5 Gyr, m_WD > 0.9 M_⊙). As a comparison, our estimate of the extra-delayed population fraction (f_extra) in this region is 6.4% (Equation (21)) and about half of them are DQs (Equation (26)). So, these number counts are consistent with the scenario that hot DQs are the evolutionary counterparts of Q-branch DQs.

In summary, based on the velocity distribution, we argue that all Q-branch DQs belong to the extra-delayed population, and these DQs account for 53 ± 16% of this population. In terms of observational properties, Q-branch DQs form a new class of DQ white dwarfs in addition to the well-understood standard-mass DQs and hot DQs. However, number counts show that hot DQs may evolve into Q-branch DQs, and both of them are likely to originate from double-WD mergers.

One additional way to test if the extra-delayed white dwarfs are also double-WD merger products is to check the wide binary fraction. The kick velocity of a few km s⁻¹ (estimated from the results of Dan et al. 2014) from a merger may destroy many wide-separation binaries, making the wide-binary fraction lower. Because the extra-delayed population is significantly enhanced on the Q branch, if the extra-delayed white dwarfs are double-WD merger products, one would expect to see a lower wide-binary fraction on the Q branch.

We cross-match the wide binaries in Gaia DR2 (El-Badry & Rix 2018, 2019) with our high-mass white dwarf sample. In the early, Q-branch, and late segments, we find 5, 4, and 7 white dwarfs with wide-separation companions out of 309, 510, and 251 white dwarfs, respectively. So, the wide-binary fraction on the Q branch is 0.8 ± 0.4%, 2σ lower than the value outside the branch (2.2 ± 0.5%). If we assume that the extra-delayed population contributes no wide-binary system, then the wide-binary fraction of the normal-cooling populations on the Q branch can be estimated as $4/[510\times (1-{F}_{\mathrm{extra}})]=(1.7\pm 0.8) \%$, consistent with the off-branch value 2.2 ± 0.5% within 1σ. Therefore, the fraction of wide binaries provides additional support for the idea that the extra-delayed white dwarfs are double-WD mergers products.