Daniel Mikkola, Paul J McMillan, David Hobbs, John Wimarsson, The velocity distribution of white dwarfs in Gaia EDR3, Monthly Notices of the Royal Astronomical Society, Volume 512, Issue 4, June 2022, Pages 6201–6216, https://doi.org/10.1093/mnras/stac434
Using a penalized maximum likelihood, we estimate, for the first time, the velocity distribution of white dwarfs in the solar neighbourhood. Our sample consists of 129 675 white dwarfs within 500 pc in Gaia Early Data Release 3. The white dwarf velocity distributions reveal a similar structure to the rest of the solar neighbourhood stars, reflecting that white dwarfs are subject to the same dynamical processes. In the velocity distribution for three magnitude-binned subsamples, we, however, find a novel structure at (U, V) = (7, −19) km s−1 in fainter samples, potentially related to the Coma Berenices stream. We also see a double-peaked feature in U − W at U ≈ −30 km s−1 and in V − W at V ≈ −20 km s−1 for fainter samples. We determine the velocity distribution and velocity moments as a function of absolute magnitude for two samples based on the bifurcation identified in Gaia Data Release 2 in the colour–magnitude diagram. The brighter, redder sequence has a larger velocity dispersion than the fainter, bluer sequence across all magnitudes. It is hard to reconcile this kinematic difference with a bifurcation caused purely by atmospheric composition, while it fits neatly with a significant age difference between the two sequences. Our results provide novel insights into the kinematic properties of white dwarfs and demonstrate the power of analytical techniques that work for the large fraction of stars that do not have measured radial velocities in the current era of large-scale astrometric surveys.
The present-day structure and the history of the Galaxy are encoded not just in the positions of its stars but also in their kinematics. It is well established that the present velocity distribution in the solar neighbourhood has a great deal of structure in it (e.g. Gaia Collaboration 2018c) for which there are multiple possible causes. Suggested origins for overdensities include dissolving open clusters, resonances from large-scale density waves such as the Galactic bar and spiral arms, accreted populations from galaxy mergers, and phase mixing from nearby satellite galaxies (e.g. Antoja et al. 2012; Kushniruk, Schirmer & Bensby 2017). Understanding this substructure is a part of understanding the dynamical history of the Milky Way.
The phase-space distribution of stars within the Milky Way has been studied extensively over the last decades (see Gaia Collaboration 2018c and references therein) to reveal this complicated structure, especially since the Hipparcos mission (Perryman et al. 1997) and more so with its successor Gaia’s (Gaia Collaboration 2016) recent second and third data release (henceforth DR2 and EDR3, respectively; Gaia Collaboration 2018a, 2021a).
The astrometry of Gaia provides proper motions and positions for ∼1.5 billion sources with great precision, which is an enormous leap forward from its predecessor, which observed ∼120 000 sources. The Gaia data have only been available for a few years but the potential for kinematic study has already been demonstrated. For example, Bovy ( 2017) accurately measured the Oort constants A and B as well as for the first time the non-axisymmetric constants C and K. In Monari et al. ( 2018), it was shown that the moving group Coma Berenices is limited to negative Galactic latitudes and likely has not undergone phase mixing in the Galactic potential. The kinematic structure of the solar neighbourhood has been studied in unprecedented detail to reveal many new and old structures (e.g, Kushniruk et al. 2017; Gaia Collaboration 2018c) as well as arches (Antoja et al. 2018). Beyond velocity space, the solar neighbourhood has also been explored in orbit space (e.g. Trick, Coronado & Rix 2019; Trick et al. 2021; Trick 2022), showing ridges which can manifest themselves as streams and structures in velocity space. Understanding the kinematic substructure of the Galaxy will require exploration in both velocity and orbit space, which becomes far more accessible due to the wealth of data provided by missions such as Gaia.
Even though EDR3 provides accurate astrometry and photometry for a great number of sources, it does not contain full 3D phase-space information for all of them as some lack measured radial velocity. This means that for most individual stars only the position and proper motions are available as in the Hipparcos catalogue. In fact, the number of sources with radial velocities in Gaia EDR3, ∼7.2 million, is dwarfed by the number of sources with at least position and proper motions, ∼1.5 billion. This means that the radial velocity sample contains only ∼0.5 per cent of sources with the full astrometric solution.
This limitation can be circumnavigated by studying properties of an entire sample rather than the individual stars. This was demonstrated in two seminal papers, Dehnen & Binney ( 1998) (hereafter DB98) and the follow-up paper Dehnen ( 1998) for Hipparcos. DB98 calculated the mean motion and velocity dispersions for populations taken from a sample of 11 865 single main-sequence stars, and Dehnen ( 1998) estimated the velocity distribution |$f(\pmb {v})$| for a similar sample of stars. Both of these papers use the approximation that the velocity distribution is consistent for the full sample and is spread across the full sky.
As previously mentioned, there are ∼1.5 billion sources in EDR3. This sample includes stars in Gaia’s G-band magnitude system as faint as G ≈ 21, which means that it contains samples of stars that have been previously unavailable for kinematic studies. This includes the a large number of white dwarfs (WDs) which, following DR2, revealed that the colour–magnitude diagram (CMD) of WDs has more structure than previously thought and displays a clear bifurcation as well as a crystallization branch (Gaia Collaboration 2018b, Tremblay et al. 2019).
Since then, there have been a few different explanations put forward and we will briefly summarize a few here. Shortly after DR2 released, El-Badry, Rix & Weisz ( 2018) used a cross-match of the 100 pc Gaia WD sample with the Montreal White Dwarf Database (MWDD; Dufour et al. 2017) to show that part of the bifurcation could be explained by an initial–final mass ratio (IFMR) that produces a bimodal WD mass distribution with peaks at both ∼0.6 and ∼0.8 M⊙ if the age distribution in multimodal. The different-mass WDs would populate different cooling tracks and produce the bifurcation. This would not produce a bifurcation in mono-age clusters, however, as massive WDs would cool before young WDs appear on the CMD.
Around the same time, Kilic et al. ( 2018) used a similar sample and showed that atmospheric composition explains the bifurcation well for 0.6 M⊙ WDs. However, they also conclude that the WD mass distribution is indeed bimodal. They suggest that this bimodality can be explained at least partly through the merger of WD binaries. More recently, Kilic et al. ( 2020) revisited the 100 pc WD sample and conducted a spectroscopic follow-up survey, thereby being able to constrain the atmospheric composition reliably. They found that the mass distribution of WDs that have H lines as the primary feature of their spectrum (DA WDs) has a sharp peak at 0.59 M⊙ with a broad shoulder, best fitted with a secondary Gaussian at 0.76 M⊙, again demonstrating the existence of a bimodal mass distribution. They test a WD model including mergers and find that it cannot produce a good fit to the observed mass distribution. They also investigate the transverse velocities of WDs, since merger products should appear as massive WDs with larger velocities than those formed from single main-sequence stars. However, the lack of young, massive WDs with large velocities coupled with the model predictions leads them to conclude that mergers are unable to explain the bimodal mass distribution of WDs. Instead, it is shown in Tremblay et al. ( 2019), Bergeron et al. ( 2019), and Kilic et al. ( 2020) that the effects of crystallization are able to create the overabundance of massive WDs in the 0.7–0.9 M⊙ range. In this scenario, the massive WDs that should already have reached the bottom of the WD sequence are subjected to cooling delays (for a detailed explanation of these effects, see e. g. Bauer et al. 2020; Blouin, Daligault & Saumon 2021).
The bifurcation in the WD CMD is well described by atmospheric differences and the bimodal WD mass distribution can arise due to core crystallization and its related effects. We explore a new direction to probe the WD bifurcated sequences using their kinematics, which could provide additional insight into the bifurcation problem. While kinematics could be affected in complex ways by the processes that produce binary merger systems, the second, fainter, sequence of the CMD should have the same velocity dispersion as the brighter sequence if it is caused by mergers or atmospheric composition. Previously, Rowell & Kilic ( 2019) performed a kinematic study of WDs in DR2 using the method of DB98 to determine the mean velocity and velocity dispersion of the WDs. However, as they only split the WD CMD in MG their results reflect a mixture of the two sequences and we expand upon this analysis by splitting the WDs across the visible bifurcation and computing the full velocity distributions in addition to velocity moments.
The paper is organized as follows: In Section 2, we briefly present the techniques of DB98 and Dehnen ( 1998) used to determine the moments and velocity distribution. The WD samples that we have used and how they are selected are presented in Section 3. Then, in Section 4 we present the moments for the bimodal sequences and in Section 5 the velocity distribution of nearby WDs is presented for the first time. The implication and significance of our results is discussed in Section 6 and our conclusion are in Section 7.
In Galactic dynamics, we decompose the space velocity into velocity towards the Galactic Centre, U, in the direction of rotation, V, and, north of the Galactic plane, W. Galactic observations, however, use a combination of the line-of-sight velocity, vr, and the combined on-sky velocity, |$\pmb {p}$| . With the release of Gaia, we have access to a large sample of stars for which positions (ℓ, b), parallax (ϖ), and proper motions (μℓ*, μb) have been measured. A subset of these will also have measurements of radial velocity vr and with these properties combined one can determine a star’s space velocity |$\pmb {v}_i$| . Without the full 3D velocity vector, we can only know a star’s tangential velocity |$\pmb {p}_i$| . Despite this, it is still possible to determine the moments of phase space, the mean velocity components, and the velocity dispersions, and we do this by making use of a deprojection technique from DB98.
This approach is very similar to that of DB98, with the only difference being that here the sample is not assumed to be perfectly isotropic. This method has seen use in, e.g., Rowell & Kilic ( 2019) for a similar sample of stars to ours.
In order to maximize the value of equation ( 25), we make use of the conjugate gradient method (e.g. Press et al. 2002). The value of the smoothing parameter α in equation ( 25) determines the balance between goodness of fit and smoothness. In order to determine an appropriate value, we use the Gaia Radial Velocity Spectrometer (RVS) sample. We create two subsamples by randomly selecting 130 000/30 000 sources (to match our all_500 and red_500 /blue_500 sample sizes) and maximize the likelihood with a range of α values to determine the optimal value given a different number of sources. Since the actual velocity distribution can be determined for the RVS sample, we can determine an appropriate α by visual inspection. We find that the values α = 10−11 and α = 3 × 10−11 appropriately reconstruct the velocity distribution and use these for our WD and red/blue samples, respectively.
We reduce the number of operations required considerably by using an approach with an increasing grid size. This method uses the initial guess of |$\pmb {\phi }$| on a crude Cartesian grid. The solution, |$\hat{\pmb {\phi }}$| , which maximizes |$\tilde{\mathscr {Q}}_\alpha$| , is then interpolated on a finer grid and used as the initial guess for a new maximization on that finer grid. We allow for up to six different grids, each grid being twice as big in each dimension. This means the initial grid, |$\pmb {n}_\mathrm{initial}$| , will be of shape |$\pmb {n}_\mathrm{final} / 2^k$| , where k is the number of grid steps and is chosen such that no dimension in |$\pmb {n}_\mathrm{initial}$| is less than 10. For our velocity distributions, we use a grid of |$\pmb {n} = [100, 100, 72]$| range. Our velocity distribution is found, following two stages of refinement, on a grid of |$\pmb {n} = [100, 100, 72]$| over the range U ∈ [−150, 150] km s−1, V ∈ [−150, 50] km s−1, and W ∈ [−80, 60] km s−1.
Our set-up is broadly consistent with the one used by Dehnen ( 1998), with a few differences. Our algorithm is built without the original rejection criterion that any star’s |$K(k|\pmb {l})$| must pass through 96 cells. This means our distributions might take into account stars that lie outside the grid. Instead of discarding these stars, we simply omit the outermost layers of the grid when we present our results to reduce the impact of numerical edge effects. Our choice of grid size and ranges also provides us with a slightly better resolution of |$\Delta \pmb {v} \sim [3, 2, 2]$| km s−1.
Since the algorithm itself does not take into account the measured uncertainties of the parameters, we statistically resample every source that we use. To do this, we draw alternative parameters for each source from a multivariate Gaussian in μα*, μδ, and ϖ, with the measured values as means and the uncertainty and correlation coefficients in the covariance matrix. We ignore the uncertainties on RA and Dec. as they are negligible. By comparing the inferred velocity distribution of these resamples with the original sample, we can estimate how significant the observed features are. We find no significant deviations from the initial sample. One example using a resampled distribution can be seen in Appendix C.
We create the WD sample by setting Mg > 9.6 + 3.7(GBP − GRP). As an additional test, we compare our WD sample to that of Gentile Fusillo et al. ( 2021) and find that 99.5 per cent of our WDs are available within their catalogue. In their paper, they calculate the probability that each source is a WD, PWD, and we find that of the cross-matched sample 98 per cent of stars have PWD > 0.9, indicating a high degree of confidence that our sources are WDs.
The WDs are split along the bimodal sequences along a line selected by eye into an upper and lower sequence (hereafter referred to as the red and blue sequences). We limit these cuts to 12 ≲ MG ≲ 14, where the bifurcation is clearest. All samples and the number of sources in them are listed in Table 1. To ensure that we are robust to the specific choice of line, we shift the line ±0.05 mag in colour, creating two different versions of the red and blue samples. We find that our results do not change when we shift the line. For WDs with d < 100 pc we show the red and blue sequences on top of the CMD in Fig. 1.
Top: CMD for the WDs within 100 pc. The colour shows the density of WDs on a 150 × 150 grid. Middle: the red and blue selections used to split across the bifurcation. The vertices of the two regions can be seen in Appendix B. Bottom: three different bins in absolute magnitude shown on top of the CMD.
Top: CMD for the WDs within 100 pc. The colour shows the density of WDs on a 150 × 150 grid. Middle: the red and blue selections used to split across the bifurcation. The vertices of the two regions can be seen in Appendix B. Bottom: three different bins in absolute magnitude shown on top of the CMD.
The names of the various samples used and the number of sources in them.
The names of the various samples used and the number of sources in them.
Our selection includes WDs at distances up to ∼500 pc. In Gaia, the nominal brightness limit is G = 20.7 (Gaia Collaboration 2021a), which in an ideal case would allow Gaia to detect sources as faint as MG = 15.7 within 100 pc and MG = 14.2 within 200 pc. In addition to this, we place a parallax uncertainty criterion of ϖ/σϖ > 10, which means that at 100 and 200 pc the uncertainties must be smaller than 1 and 0.5 mas, respectively. In Gaia, the typical parallax uncertainty at G = 20.7 for five-parameter solutions is 1.3 mas (Gaia Collaboration 2021a). The brightness and uncertainty limits mean that beyond 100 pc we will start to be affected by incompleteness, and therefore Malmquist bias. Conversely, our samples within 100 pc are free of this bias, especially the red and blue samples, which only go as faint as MG = 14. We also use two WD samples limited in distance to 200 and 500 pc, which are limited by this bias. The effect of the bias will not be exactly the same for the two sequences as the blue sequence is slightly fainter than the red. However, the shift in magnitude between the sequences is sufficiently small that we can neglect it and assume they share the same bias. This means that any difference we see between the two sequences is not be due to the Malmquist bias as it affects the sequences in effectively the same way. A two-sample Kolmogorov–Smirnov (KS) test on the velocity dispersion of the samples between 0 and 100 pc and between 100 and 200 pc shows that these are drawn from the same population, so we can be confident that Malmquist bias has not had a significant effect on our sample out to 200 pc (see the following).
We also split the sample into three bins in absolute magnitude. These bins are chosen such that they have an equal number of WDs from the all_500 sample. The bins are illustrated in Fig. 1, labelled A , B , and C , and represent a brighter, intermediate, and faint selection of WDs, respectively. As the WDs grow older, they will cool and become fainter, so by looking at different magnitudes we should be able to detect age-dependent variation in the velocity distribution.
In Section 2.1, we showed how we can calculate the velocity dispersion of a sample of stars using only their positions and tangential velocities. We calculate 3D velocity dispersions for a moving window in absolute magnitude separately for the WD samples limited to 100 and 200 pc. We estimate the uncertainty within the windowed sample by calculating the moments for 500 bootstrapped samples. The standard deviation of all of the resulting statistics is used to provide a 1σ uncertainty. The complete 100 pc sample’s velocity dispersions are seen in Fig. 2a. Generally, we would expect the velocity dispersion to increase with fainter magnitude due to dynamical heating over time. This is true overall, though the component σV is approximately constant for the red and black sequences over the range of magnitudes considered.
(a) Dispersions in U, V, and W calculated for samples all_100 , red_100 , and blue_100 using a moving window in absolute magnitude and shown with black, red, and blue colours, respectively. The shaded region shows the 1σ uncertainty. The red and blue lines appear to be separate for brighter WDs in all three directions but become mixed towards the fainter end of the sequence. (b) Same as (a) but for the samples all_200 , red_200 , and blue_200 . For these WDs, the split between the red and the blue sequences is much more pronounced and now clearly so at all absolute magnitudes. The red sequence appears to have a larger velocity dispersion in all directions and at almost all absolute magnitudes.
(a) Dispersions in U, V, and W calculated for samples all_100 , red_100 , and blue_100 using a moving window in absolute magnitude and shown with black, red, and blue colours, respectively. The shaded region shows the 1σ uncertainty. The red and blue lines appear to be separate for brighter WDs in all three directions but become mixed towards the fainter end of the sequence. (b) Same as (a) but for the samples all_200 , red_200 , and blue_200 . For these WDs, the split between the red and the blue sequences is much more pronounced and now clearly so at all absolute magnitudes. The red sequence appears to have a larger velocity dispersion in all directions and at almost all absolute magnitudes.
Cumulative distribution of the Q statistic described in the text for the red and blue sequences when comparing samples limited to within 100 pc or between 100 and 200 pc. The black line shows the Gaussian CDF with μ = 0 and σ = 1.
Cumulative distribution of the Q statistic described in the text for the red and blue sequences when comparing samples limited to within 100 pc or between 100 and 200 pc. The black line shows the Gaussian CDF with μ = 0 and σ = 1.
If the bifurcation is caused by a bimodal WD mass distribution, this could be explained by the progenitors of the WDs in the two different sequences. The lower mass WDs will come from less massive main-sequence stars, which would mean that these progenitors will have been dynamically heated for a longer duration of time and the lower mass WDs would be born with larger dispersions.
The separation between red and blue sequences is larger at brighter magnitudes, which correspond to younger WDs. This can be expected if the blue sequence comes from massive WDs that are younger than their red-sequence, less massive, counterparts. Most of the dynamical heating occurs during the first few Gyr of a stars life (e.g. Nordström et al. 2004; Binney & Tremaine 2008) and so would already have occurred for the red-sequence WDs but still be occurring for the blue-sequence ones. Over time as WDs on both sequences cool to fainter magnitudes, their heating tracks would align to become more parallel, which is visible in our figures.
Undeniably, the bifurcation in the CMD is well described by a difference in the atmospheric composition of WDs (Bédard et al. 2021). It is however difficult to reconcile with the observed difference in the kinematics that we display here, which requires significantly different dynamical heating histories between the two sequences. If there is a significant age difference between the two populations, however, this would result from the expected heating from the Galactic disc. We can estimate the age difference of stars in the two sequences using the age–velocity dispersion relationship from Aumer & Binney ( 2009) and applying it to the total velocity dispersion of the two samples, red_500 and blue_500 . This suggests a typical age of ∼7.7 Gyr and ∼4.9 Gyr respectively for the red and blue sequences. We are unaware of any mechanism that could explain these kinematics purely by a difference in atmospheric composition.
As mentioned in Section 3, we also selected WDs out to ϖ = 2 mas or 500 pc. Beyond 200 pc, the separation between red and blue sequences remains the same and did not provide any further insights.
We apply the maximum penalized likelihood estimate to all of the samples to produce a full 3D velocity distribution for which we show the projections in the planes U–V, U–W, and V–W. To identify where the peaks of the distribution are, we apply the peak_local_max function available as part of the scikit-image 1python package (van der Walt et al. 2014). The results in the U–V plane are presented first in Fig. 4. In the first row, we compare the velocity distribution of all our WDs with those of WDs identified as belonging to the red and blue sequences as laid out in Section 3. It is immediately clear that the WD velocity distribution is overall very similar to the rest of the solar neighbourhood (e.g. Gaia Collaboration 2018c), which is not too surprising considering that the WDs are older stars of the same population.
The velocity distribution of WDs in U and V. The top row shows the distribution for the whole WD sequence and the red and blue sequences corresponding to samples all_500 , red_500 , and blue_500 in Table 1. Contour lines are constructed so as to contain 95, 90, 80, 68, 50, 33, 21, 12, 6, and 2 per cent of all sources in the sample and magenta crosses show identified peaks with a peak-finding algorithm described in Section 5. The second row contains samples A, B, and C , as indicated, which are magnitude bins in the full sample. Three horizontal arches are illustrated in sample C and explained in Section 5. The third row takes the distributions for the three samples A, B, and C and subtracts the mean of all three samples. The contour lines are the same as the individual distributions. The fourth row shows the first nine groups identified in Antoja et al. ( 2012) as black crosses for comparison.
The velocity distribution of WDs in U and V. The top row shows the distribution for the whole WD sequence and the red and blue sequences corresponding to samples all_500 , red_500 , and blue_500 in Table 1. Contour lines are constructed so as to contain 95, 90, 80, 68, 50, 33, 21, 12, 6, and 2 per cent of all sources in the sample and magenta crosses show identified peaks with a peak-finding algorithm described in Section 5. The second row contains samples A, B, and C , as indicated, which are magnitude bins in the full sample. Three horizontal arches are illustrated in sample C and explained in Section 5. The third row takes the distributions for the three samples A, B, and C and subtracts the mean of all three samples. The contour lines are the same as the individual distributions. The fourth row shows the first nine groups identified in Antoja et al. ( 2012) as black crosses for comparison.
We can identify a few familiar moving groups such as the Hyades, Pleiades, and Hercules across all three samples very clearly. To a lesser extent, we can also identify higher velocity moving groups beyond U > 25 km s−1 such as Wolf 630 and Dehnen98.The moving groups in the blue sequence are more narrowly distributed around their group mean overall than those of the red sequence, where the groups are spread out more across the central regions in the space due to the overall larger velocity dispersion shown in Section 4.
We show the three samples of different absolute magnitudes in the second row of Fig. 4, going from brightest (A ) to faintest (B ). For the individual distributions, we see the same type of structure as we did for the whole sample. It appears as though there is more structure in these subsamples than in the former three. However, as the samples contain less stars than the all_100 sample but use the same smoothing parameter α they have a noisier distribution. The fainter sample, C , has older stars and thus a larger velocity dispersion, as expected, which ‘smears’ out the distribution, which reveals an arch-like structure in C , with three horizontal arches at V of 0, −20, and −40 km s−1, which are illustrated in the plot. This type of structure has been shown to exist in the solar neighbourhood (Gaia Collaboration 2018c) and can be attributed to dynamical resonances with the spiral arms and Galactic bar (Trick et al. 2019).
The third row takes the three samples A , B , and C , and subtracts the mean of all three to emphasize the parts of the distribution that are shared between the samples or are independently significant. Our expectation is that the bright, young sample with low velocity dispersion (A ) should be overabundant near the centre of the distribution and the faint sample should show the opposite. Stars with very large |V| belong to populations with their guiding centres at significantly larger/smaller Galactocentric radii and would be visiting the solar neighbourhood when their orbits are sufficiently dynamically heated. Thus, they would be older and part of a fainter sample. This effect can be seen in the third row, where samples A and B dominate the centre of the distribution while C is represented more in the wings with, for example, the Hercules moving groups being overabundant. Asymmetric drift also causes the mode of the distribution to shift towards more negative V.
The region around (U, V) ≈ (7, −19) km s−1 is underdense in A and B but surprisingly is overabundant in C , with an identified peak in that region as well. A review of recent articles that investigate the velocity distribution in U − V for the solar neighbourhood reveals that Antoja et al. ( 2012) and Gaia Collaboration ( 2018c) do not identify a known moving group in this region, while Kushniruk et al. ( 2017) does and attributes it to the Coma Berenices stream.
The distributions in U − W and V − W in Figs 5 and 6 are not as rich in structure as the U− V plane, but we can, however, identify Coma Berenices and Pleiades by comparing the two planes. Our choice of smoothing parameter α is chosen to fit best with the sample all_500 and hence will be underestimated for the smaller samples (red_500 and blue_500 ) and as such the features seen in them may disappear with an α that is appropriate for smaller samples. A feature identified by Dehnen ( 1998) is a double-peaked feature along W in the U − W plane at U ≈ −30 km s−1 and in the V− W plane at V ≈ −20 km s−1. This double-peaked feature is very vaguely present in U− W for sample C . In the V− W plane, however, we can identify the feature clearly in both red_500 and blue_500 , and perhaps vaguely in B and C . We mark the proposed double-peak feature with circles around the involved features. This seems to imply that the double-peaked feature, of which the Pleiades appears a part, is limited to fainter, possibly older stars.
Same as Fig. 4 but in U and W. The circles in panel C mark the peaks that are suggested to be part of a double-peaked feature discussed in Section 5.
Same as Fig. 4 but in U and W. The circles in panel C mark the peaks that are suggested to be part of a double-peaked feature discussed in Section 5.
Same as Fig. 5 but in V and W.
Same as Fig. 5 but in V and W.
In Table 2, we summarize the features we can see and mark whether we can see them clearly, somewhat, or not at all. We also list their locations in all three velocity dimensions (see for comparison table 2 from Dehnen 1998). For most of the moving groups, we are able to associate a previously known one following Kushniruk et al. ( 2017), but no known group accurately describes feature 5 in the table, meaning that it may be a new feature. Compared with previous results, the features presented here are located mostly at lower velocities (|v| < 50 km s−1). It may be that the WDs in groups at larger velocities are not numerous enough to appear in the distributions or that their dispersion is too great.
Identified features in the velocity distributions at approximate coordinates. We omit the velocity in W when it cannot be determined accurately. Filled circles and hollow circles are for clear and weak features, respectively. Group 5 has an asterisk indicating unclear membership in the moving group, which is further discussed in Section 6.2.
Identified features in the velocity distributions at approximate coordinates. We omit the velocity in W when it cannot be determined accurately. Filled circles and hollow circles are for clear and weak features, respectively. Group 5 has an asterisk indicating unclear membership in the moving group, which is further discussed in Section 6.2.
In his analysis of Hipparcos data, Dehnen ( 1998) found evidence for the stellar warp by investigating 〈W〉 as a function of V. For V > 10 km s−1, the velocity distribution was skewed towards positive 〈W〉. We perform the same analysis here of our largest WD sample all_500 and the results can be seen in Fig. 7. To do this, we require a larger grid of |$\pmb {n} = [200, 152, 152]$| covering U ∈ [−200, 200] km s−1, V ∈ [−150, 150] km s−1, and W ∈ [−150, 150] km s−1 to ensure that none of the edge effects discussed in Section 2.2.1 affect the results. We find that our results also show increasing 〈W〉 with increasing V.
Mean of the vertical motion, 〈W〉, in the central regions of f(V, W) for sample all_500 . The location of the LSR is marked by a triangle and is taken from Schönrich, Binney & Dehnen ( 2010). The error bars show the 1σ Poisson noise. The dashed black line shows a weighted linear fit.
Mean of the vertical motion, 〈W〉, in the central regions of f(V, W) for sample all_500 . The location of the LSR is marked by a triangle and is taken from Schönrich, Binney & Dehnen ( 2010). The error bars show the 1σ Poisson noise. The dashed black line shows a weighted linear fit.
We compare our results to Schönrich & Dehnen ( 2018), who use the Tycho-Gaia Astrometric Solution (TGAS) to determine the average vertical velocity as a function of azimuthal velocity and angular momentum for stars in the two cones towards the Galactic Centre and anti-Centre with angular radius of 30° (their fig. 6). Since our samples are in the solar neighbourhood, the velocity in V is a proxy for angular momentum since Lz = −R(V⊙ + V). Just as in our sample, the average vertical velocity increases with angular momentum and shows a dip just beyond the local standard of rest (LSR). At large angular momentum, corresponding to V > 15 km s−1 in the solar neighbourhood, there is a significant increase in 〈W〉 for all of their samples as well as ours. Towards V = 50 km s|$^{-1}\, \langle W\rangle$| decreases again, which is consistent with a similar decrease at higher angular momentum seen by Schönrich & Dehnen ( 2018).
More recently, the Gaia Collaboration ( 2021b) also looked at the vertical velocity profile of stars outside the solar radius as a function of angular momentum. We see this in their fig. 11, which shows this profile for several different stellar types, of which the young main-sequence stars are the ones to show a decrease in vertical velocity at |$L_z^\ast \sim 2400$| kpc km s−1 corresponding to around V = 60 km s−1. However, our sample and the Gaia TGAS sample are in the solar neighbourhood, where stars with such large angular momentum must have large radial excursions to be included in the sample, and therefore are very likely to be old. This is in contrast to the stars in the Gaia Collaboration ( 2021b), which are young and located beyond 10 kpc from the Galactic Centre. While the decrease in the vertical velocity profile we see could therefore be related to these decreases seen in other data sets, we note that the decrease is still within 1σ uncertainty, and therefore, the most we can say is that the results are mutually consistent.
For V < 20 km s−1, where the Poisson noise is not as great, the average velocity, 〈W〉, has a relatively small tilt. We perform a weighted linear fit for V ∈ [−80, 50] km s−1 and find a slope of ∼0.03 (shown in Fig. 7), which is somewhat larger than the ∼0.02 found by Schönrich & Dehnen ( 2018). If we instead restrict our fit to the relatively well-constrained region V ∈ [−80, 20] km s−1, we find a slope of ∼0.025, in agreement with the previous results.
In Section 1, we reviewed the literature on the nature of the bifurcation of the CMD in the Gaia data. We know that there exists a bimodality in the mass distribution of spectral class DA WDs (Kilic et al. 2018, 2020; Jiménez-Esteban et al. 2018) at around ∼0.6 M⊙ and ∼0.8 M⊙.
Cooling tracks for these masses of DA WDs would produce a bifurcation that appears very similar to the one visible in Gaia (El-Badry et al. 2018). However, massive WDs will form sooner and begin cooling at an earlier time than the less massive ones, which means that for these cooling tracks to be simultaneously visible in the CMD the massive WDs must be (1) formed through the merger of lower mass WD binaries, (2) formed at a later time than the less massive WDs, or (3) cooling slower due to crystallization. Scenario (1) was suggested by Kilic et al. ( 2018) but subsequently dismissed by Kilic et al. ( 2020) due to merger models being unable to produce a mass distribution that fit the observations as well as, and, perhaps more importantly, not being able to find a significant amount of young and massive DA WDs with high velocities. This leaves scenarios (2) and (3). In the first of these scenarios, a multimodal age distribution could produce massive WDs at a later time for a bimodal mass distribution that would have a smaller velocity dispersion as they have not had as much time to be heated. In the latter scenario, crystallization (Tremblay et al. 2019) will result in a slowdown of the cooling of massive WDs, which means the massive sequence can consist of a mixture of young and older massive WDs. Both of these scenarios would cause the massive WDs to have a lower velocity dispersion. In Fig. 2, it is clear that this behaviour is observed. The velocity distributions for the red and blue sequences in Figs 4, 5, and 6 exhibit this as well, with the red sequence showing a larger velocity dispersion.
Atmospheric composition can elegantly explain the visible bifurcation in the Gaia data with an upper DA branch and a lower branch of He-dominated WDs with trace amounts of H or other metals. The bimodal WD mass distribution fits well with this explanation when described by the process of crystallization or by the inclusion of young massive WDs. However, atmospheric composition alone does not provide the difference in kinematics between the two sequences that we identify. This difference arises naturally, however, with a bimodal WD mass distribution from a multimodal age distribution.
The cooling tracks of ∼0.6 M⊙ and ∼0.8 M⊙ DA WDs would also lie in the region where the bifurcation exists (e.g. Bergeron et al. 2019). We cross-match our red_100 and blue_100 with the MWDD to determine what fractions of stars are DA or non-DA and find that around 85 per cent of the red sample cross-match and 39 per cent of the blue sample cross-match are DAs. The red cross-match shows that the this sequence is likely to be comprised of older, less massive DAs. The DAs in the blue cross-match have higher mass with a mean of ∼0.8 M⊙ as expected in keeping with previous results (e. g, Kilic et al. 2020).
The bifurcation seen in Gaia clearly has contributions from both atmospheric differences and different WD mass cooling tracks. It remains unclear whether the origin of the massive WDs is recent bursts of star formation or a pile-up of WDs with a mixture of ages due to delayed cooling from crystallization, or a mixture of both. Fantin et al. ( 2019) investigate the star formation history of the Galactic disc using WDs and suggest that star formation increases 3.3 ± 1.8 Gyr ago and is roughly constant for ∼5 Gyr prior to that. A detailed study of the age distribution of the WD population using accurately determined ages, masses, and spectral types would provide valuable insight into these questions and analytical modelling of WD formation and evolution following bursts of star formation would be a good avenue to test the formation avenues of these WDs.
For the first time, we present the velocity distribution of WDs in the solar neighbourhood in addition to their velocity moments for such a large sample. We find that the WD velocity distribution in (U, V) shares many features with the velocity distribution of main-sequence stars when comparing our results with recent maps of the kinematic structure of the solar neighbourhood like those of the Gaia Collaboration ( 2018c), Kushniruk et al. ( 2017), or Antoja et al. ( 2012). When we separate our sample along this bifurcation, we can see very similar velocity distributions, with the notable differences being the red sample having a larger velocity dispersion.
The CMD is also split into three equally sized samples based on their absolute magnitude. The mean velocity distribution of all three subsamples is subtracted from each individual subsample, which reveals an unexpected overdensity for the faintest sample, C , in the region (U, V) ≈ (7, −19) km s−1. The location does not match conclusively to any of the known moving groups and is only identified in Kushniruk et al. ( 2017), who attribute it to be a part of the Coma Berenices moving group along with two other identified groups. If this feature is limited to fainter, older stars, it could suggest its origin is dynamical. In Monari et al. ( 2018), it is shown that the Coma Berenices moving group is not vertically phase mixed and is localized to negative b only, and this is suggested to be due to a recent passing by a dwarf galaxy such as Sagittarius, which fits well with passages suggested in the literature (for a summary, see the lower panel of fig. 2 of Ruiz-Lara et al. 2020).
The double-peaked feature identified in W is limited to fainter stars, which should be part of an older sample. The feature is symmetrical around the average vertical velocity, which suggests that this feature might be in the brighter sample as a single feature around the mode containing younger, less dynamically heated stars.
Beyond smaller scale structure, we find that the velocity distribution is very similar to that of main-sequence stars, which can be expected since WDs are subject to the same dynamical processes as other stars. The distributions also show very clear arches like those seen in dynamical studies of main-sequence stars (e.g. Trick et al. 2019). To further narrow down the origin of the observed features requires comparison with model predictions or investigating the ages and abundances of their associated stars. For example, It would be possible to investigate the age evolution of dynamical features in the distributions by studying their appearance for several absolute magnitude bins and associating the absolute magnitude with ages using WD cooling tracks.
We present the velocity distributions of WDs in the solar neighbourhood, which, to our knowledge, has not been done previously. The velocity distributions are estimated through a penalized maximum likelihood. The data we use come from Gaia EDR3 and are filtered to select a clean and unbiased sample of WDs. We split the WD CMD across the established bifurcation and into several absolute magnitude bins and find that the velocity distribution is similar to that of main-sequence stars from previous studies and displays many well-known moving groups. We also identify a novel structure located at (U, V) ≈ (7, −19) km s−1, which appears only in bins that are fainter or possibly older, and discuss possible explanations for the feature. We also identify a double-peaked feature in W, previously established in the Hipparcos data by Dehnen ( 1998), involving mostly fainter stars.
We also explore the mean velocities and velocity dispersions as a function of absolute magnitude and compare them between the two established sequences in the Gaia WD CMD. We find that the brighter sequence has larger velocity dispersion than the faint one across all magnitudes: The sequences are two separate kinematic populations. This cannot be explained if the origin of the second sequence is due to WD binary mergers or solely atmospheric composition. Our results are consistent with the observed WD bimodal mass distribution with a multimodal age distribution.
The results of our study shed light on the bifurcation in the Gaia WD CMD and explore the possibility of accessing the majority of sources in Gaia that lack radial velocities. We plan to use EDR3 to investigate the velocity distribution of the solar neighbourhood using as many stars as we can include to probe the kinematic structure of the Milky Way for new insights.
This work 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. This research made use of astropy , 2 a community-developed core python package for astronomy (Astropy Collaboration 2013, 2018).
We thank members of the Lund Observatory for helpful comments and ideas. We give special thanks to Ross Church for constructive discussions and comments. Computations for this study were performed on equipment funded by a grant from the Royal Physiographic Society in Lund. PM was supported by project grants from the Swedish Research Council (vetenskaprådet reg: 2017-03721 and 2021-04153). DH and PM gratefully acknowledge support from the Swedish National Space Agency (SNSA Dnr 74/14 and SNSA Dnr 64/17).
All the data analysed in this paper are publicly available from the Gaia EDR3 archive (http://gea.esac.esa.int/archive/). The 3D probability distributions used in Figs 4, 5, and 6 are available upon request from the corresponding author.
https://gea.esac.esa.int/archive/
Antoja T. et al. , 2012, MNRAS, 426, L1 10.1111/j.1745-3933.2012.01310.x
Antoja T. et al. , 2018, Nature, 561, 360 10.1038/s41586-018-0510-7
Astropy Collaboration , 2013, A&A, 558, A33 10.1051/0004-6361/201322068
Astropy Collaboration , 2018, AJ, 156, 123 10.3847/1538-3881/aabc4f
Aumer M. , Binney J. J. , 2009, MNRAS, 397, 1286 10.1111/j.1365-2966.2009.15053.x
Bauer E. B. , Schwab J. , Bildsten L. , Cheng S. , 2020, ApJ, 902, 93 10.3847/1538-4357/abb5a5
Bédard A. , Brassard P. , Bergeron P. , Blouin S. , 2021, preprint (arXiv:2112.09989)
Bergeron P. , Dufour P. , Fontaine G. , Coutu S. , Blouin S. , Genest-Beaulieu C. , Bédard A. , Rolland B. , 2019, ApJ, 876, 67 10.3847/1538-4357/ab153a
Binney J. , Tremaine S. , 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton
Blouin S. , Daligault J. , Saumon D. , 2021, ApJ, 911, L5 10.3847/2041-8213/abf14b
Bovy J. , 2017, MNRAS, 468, L63 10.1093/mnrasl/slx027
Dehnen W. , Binney J. J. , 1998, MNRAS, 298, 387 (DB98) 10.1046/j.1365-8711.1998.01600.x
Dufour P. , Blouin S. , Coutu S. , Fortin-Archambault M. , Thibeault C. , Bergeron P. , Fontaine G. , 2017, in Tremblay P. E. , Gaensicke B. , Marsh T. , eds, ASP Conf. Ser. Vol. 509, 20th European White Dwarf Workshop. Astron. Soc. Pac, San Francisco, p. 3
El-Badry K. , Rix H.-W. , Weisz D. R. , 2018, ApJ, 860, L17 10.3847/2041-8213/aaca9c
Fantin N. J. et al. , 2019, ApJ, 887, 148 10.3847/1538-4357/ab5521
Gaia Collaboration , 2016, A&A, 595, A1 10.1051/0004-6361/201629272
Gaia Collaboration , 2018a, A&A, 616, A1 10.1051/0004-6361/201833051
Gaia Collaboration , 2018b, A&A, 616, A10 10.1051/0004-6361/201832843
Gaia Collaboration , 2018c, A&A, 616, A11 10.1051/0004-6361/201832865
Gaia Collaboration , 2021a, A&A, 649, A1 10.1051/0004-6361/202039657
Gaia Collaboration , 2021b, A&A, 649, A8 10.1051/0004-6361/202039714
Gentile Fusillo N. P. et al. , 2021, MNRAS, 508, 3877 10.1093/mnras/stab2672
Jiménez-Esteban F. M. , Torres S. , Rebassa-Mansergas A. , Skorobogatov G. , Solano E. , Cantero C. , Rodrigo C. , 2018, MNRAS, 480, 4505 10.1093/mnras/sty2120
Kilic M. , Hambly N. C. , Bergeron P. , Genest-Beaulieu C. , Rowell N. , 2018, MNRAS, 479, L113 10.1093/mnrasl/sly110
Kilic M. , Bergeron P. , Kosakowski A. , Brown W. R. , Agüeros M. A. , Blouin S. , 2020, ApJ, 898, 84 10.3847/1538-4357/ab9b8d
Kushniruk I. , Schirmer T. , Bensby T. , 2017, A&A, 608, A73 10.1051/0004-6361/201731147
Lindegren L. , 2018, Re-normalising the astrometric chi-square in Gaia DR2, GAIA-C3-TN-LU-LL-124. http://www.rssd.esa.int/doc_fetch.php?id = 3757412
Monari G. et al. , 2018, RNAAS, 2, 32 10.3847/2515-5172/aac38e
Nordström B. et al. , 2004, A&A, 418, 989 10.1051/0004-6361:20035959
Perryman M. A. C. et al. , 1997, A&A, 500, 501
Press W. H. , Teukolsky S. A. , Vetterling W. T. , Flannery B. P. , 2002, Numerical recipes in C+ + : the art of scientific computing,Cambridge University Press, Cambridge, UK
Riello M. et al. , 2021, A&A, 649, A3
Rowell N. , Kilic M. , 2019, MNRAS, 484, 3544 10.1093/mnras/stz184
Ruiz-Lara T. , Gallart C. , Bernard E. J. , Cassisi S. , 2020, Nat. Astron., 4, 965 10.1038/s41550-020-1097-0
Schönrich R. , Dehnen W. , 2018, MNRAS, 478, 3809 10.1093/mnras/sty1256
Schönrich R. , Binney J. , Dehnen W. , 2010, MNRAS, 403, 1829 10.1111/j.1365-2966.2010.16253.x
Tremblay P.-E. et al. , 2019, Nature, 565, 202 10.1038/s41586-018-0791-x
Trick W. H. , Coronado J. , Rix H.-W. , 2019, MNRAS, 484, 3291 10.1093/mnras/stz209
Trick W. H. , Fragkoudi F. , Hunt J. A. S. , Mackereth J. T. , White S. D. M. , 2021, MNRAS, 500, 2645 10.1093/mnras/staa3317
van der Walt S. et al. , 2014, PeerJ, 2, e453
The following query has been used on the Gaia archive 3 and is detailed in Section 3:
select bp_rp, phot_g_mean_mag, phot_bp_rp_excess_factor, ra, dec, parallax, pmra, pmdec ,
phot_bp_rp_excess_factor - (1.162004 + 0.011464* bp_rp
phot_bp_rp_excess_factor - (1.154360 + 0.033772* bp_rp
phot_bp_rp_excess_factor - (1.057572 + 0.140537*bp_rp)
phot_bp_rp_excess_factor
and parallax_over_error > 10
Table B1 lists the vertices for the intersect between red and blue regions of the WD CMD outlined in Section 3.
Vertices for the regions in MG and GBP − GRP that make up our red and blue WD sequences.
Vertices for the regions in MG and GBP − GRP that make up our red and blue WD sequences.
At the end of Section 2.2, we explained how we used a statistical resample of our proper motions and parallaxes with their measured uncertainties to estimate the effect that the uncertainties might have. As they are not significant, we omit them from the main paper and show a single example of one of these resamples in Figs C1– C3.
Same as Fig. 4 using a statistical resample of the proper motions and parallaxes.
Same as Fig. 4 using a statistical resample of the proper motions and parallaxes.
Same as Fig. 5 using a statistical resample of the proper motions and parallaxes.
Same as Fig. 5 using a statistical resample of the proper motions and parallaxes.
Same as Fig. 6 using a statistical resample of the proper motions and parallaxes.
Same as Fig. 6 using a statistical resample of the proper motions and parallaxes.
Oxford University Press is a department of the University of Oxford. It furthers the University's objective of excellence in research, scholarship, and education by publishing worldwide
Sign In or Create an Account
This PDF is available to Subscribers Only
For full access to this pdf, sign in to an existing account, or purchase an annual subscription.