Inferring properties of the local white dwarf population in astrometric and photometric surveys
The Gaia mission will provide precise astrometry for an unprecedented number of white dwarfs (WDs), encoding information on stellar evolution, Type Ia supernovae progenitor scenarios, and the star formation and dynamical history of the Milky Way. With such a large data set, it is possible to infer properties of the WD population using only astrometric and photometric information. We demonstrate a framework to accomplish this using a mock data set with SDSS ugriz photometry and Gaia astrometric information. Our technique utilises a Bayesian hierarchical model for inferring properties of a WD population while also taking into account all observational errors of individual objects, as well as selection and incompleteness effects. We demonstrate that photometry alone can constrain the WD population’s distributions of temperature, surface gravity and phenomenological type, and that astrometric information significantly improves determination of the WD surface gravity distribution. We also discuss the possibility of identifying unresolved binary WDs using only photometric and astrometric information.
keywords:white dwarfs – stars: statistics – astrometry
White dwarfs (WDs) are the remnants of stars with initial masses M (Ritossa et al., 1996; Smartt et al., 2009). The local WD population carries information about the Galaxy’s star formation and dynamical history, and constrains models of stellar evolution (Winget et al., 1987; García-Berro & Oswalt, 2016; El-Badry et al., 2018).
The Sloan Digital Sky Survey (SDSS) catalogued spectroscopically confirmed WDs (Kleinman et al., 2013; Kepler et al., 2015). A fundamental difficulty in studying WDs is that their mass is degenerate with distance. This degeneracy can be broken with high quality spectrometry and accurate atmospheric models. The Gaia mission, which recently published its second data release (DR2), is expected to increase the number of known WDs by approximately an order of magnitude (Jordan, 2007; Carrasco et al., 2014); Kruijssen et al. (2018) and Gentile Fusillo et al. (2018) have recently published WD catalogues, the latter containing 260,000 high-confidence WDs. Gaia also provides astrometric information for local neighborhood WDs. For comparison, the astrometric mission Hipparcos had a limiting apparent magnitude of (Perryman et al., 1997), while Gaia will see objects as faint as (with this limit, a WD with a mass of 0.6 M and effective temperature of 8,000 K is seen to pc, assuming no dust extinction).
In a model of the WD population, it is physically meaningful to divide the total population into WD sub-populations. WDs form a family of atmospheric types, where the main classification is between DA and DB, depending on if the envelope is hydrogen- or helium-dominated (Tremblay & Bergeron, 2008; Bergeron et al., 2011; Koester & Kepler, 2015). DA and DB stars can be identified with accurate photometry, as demonstrated by Mortlock et al. (2009). The halo WD population is kinematically warmer and older than the disk WD population, such that inferring and comparing properties of these sub-populations can yield information on the star formation and kinematic history of our Galaxy (Isern et al., 1998; Dame et al., 2016). The sub-population of binary WD systems holds information about stellar evolution (Postnov & Yungelson, 2014) and Type Ia supernovae progenitor scenarios (Livio & Mazzali, 2018), but unresolved binaries are very difficult to identify even with high-quality spectroscopy.
With the enormous size of the Gaia data set, there is great potential for inferring properties of the WD population using photometry and astrometry, rather than the smaller data set of spectroscopically observed WDs. In this work we demonstrate how to infer properties of the WD population in the solar neighborhood, using SDSS ugriz photometry and Gaia astrometry. We generate a mock data sample of WDs from a population model of temperature, surface gravity, and spatial number density distribution, of DA and DB atmospheric types. All sample objects have photometry and parallax information with observational errors expected from SDSS and Gaia, and sample construction selection effects are taken into account. We also discuss the possibility of identifying binary WD systems and demonstrate how to do so using photometric and astrometric information alone.
This paper is organized as follows. We outline our model for the WD population and the observational data that we consider, in Sec. 2 and Sec. 3 respectively. We present out method of statistical inference in Sec. 4. followed by Sec. 5, where we generate a mock data catalogue and infer the model parameters from that data. We discuss possible extensions to the WD model in Sec. 6, such as differentiating between disk and halo sub-populations, as well as the possibility of identifying unresolved double-degenerate binary WD systems. Finally, in Sec. 7 we present our conclusions.
We model the Milky Way’s population of WDs by a spatial distribution, and distributions in intrinsic WD properties. A WD is parametrized by effective temperature (), surface gravity (), phenomenological type (), and spatial position (). These are listed in Table 1, along with the population parameters and data.
There are five population parameters in our model, encapsulated in a vector : the population parameter , which parametrizing the distribution of temperatures; , , which parametrize the distribution of surface gravities; and which is the fraction of DB-type WDs.
The distribution of effective temperature is parametrized as
where is the Heaviside step function, and K and ,000 K the lower and upper bounds to the effective temperature.
The distribution of surface gravity is parametrized by
where and are the lower and upper bounds to the surface gravity, and is Student’s -distribution of width and variance , such that the quantity characterizes the heaviness of the distribution’s tails. 111 Student’s -distribution is defined (3) where we reparametrize the shape parameter according to (4)
The type of the object constitutes a third parameter of the intrinsic WD properties, called . This is a label, denoting for example if the WD is of DA or DB atmospheric classification. The probabilities are written in terms of the fraction of DB stars, as
Given the intrinsic properties of a WD, the absolute magnitude in different photometric bands can be calculated using a stellar model. In this paper, we use the Bergeron atmospheric models for WDs (Bergeron et al., 1995; Finley et al., 1997; Bergeron et al., 2001; Fontaine et al., 2001).
Also included in our model is a WD number density function, based on a Galactic model by Jurić et al. (2008), expressed in terms of cylindrical coordinates , and :
where , , , , , , , and . Assuming a solar position of (where the height of the Sun above the plane is neglected), the Galactic coordinates are given by cylindrical heliocentric coordinates through
The azimuthal angle can be neglected, as the Galaxy is assumed to be axisymmetric in this model.
|slope of temperature distribution|
|mean surface gravity|
|width of distribution|
|heaviness of the tails of the distribution|
|the fraction of DB type WDs|
|phenomenological type (DA or DB)|
|observed photometric magnitude|
|observed Galactic longitude|
|observed Galactic latitude|
The data for a given WD are apparent magnitude measurements in photometric bands (), a parallax measurement (), angular position ( and ), including the error models of these observables.222A hatted quantity (e.g. ) refers to an observed value, while a non-hatted quantity (e.g. ) refers to an observable’s true value. The angles and are written without hats, as their measurement uncertainties are so small that they can be neglected. The data characterizing a WD with index , called , is listed in Table 1.
The likelihood associated with a WD is
where is the normal distribution of mean and standard deviation . The factor containing parallax information is dropped when no parallax information is present. The apparent magnitudes are functions of the object parameters, coming from a stellar model.
In this work, we use SDSS photometry in ugriz colour bands, and a Bergeron atmospheric stellar model, as discussed in Sec. 2. In order to assign realistic uncertainties to the mock data that we generate, we use median uncertainties of the SDSS DR9 catalogue, in bins of observed apparent magnitude of width 0.5 mag. These median uncertainties are shown in Fig. 1. We limit the minimum ugriz magnitude uncertainty to 0.01 mag, in order to account for possible systematic uncertainties in the Bergeron atmospheric model.
We use parallax information with the precision of Gaia DR2, which has a limiting magnitude of . As listed in Lindegren et al. (2018), the parallax uncertainties of Gaia DR2 are about mas for sources with , about mas for sources with , and about mas at . In order to account for this magnitude dependence, we interpolate these points as shown in Fig. 2. The errors in the Gaia photometric -band are interpolated in the same manner as for the parallax, with errors of 0.3 mmag for , 2 mmag for , and 10 mmag for .
4 Statistical model
A directed acyclical graph (DAG) of our statistical model is shown in Fig. 3. In the DAG, the constants and parameters of the model are inscribed in squares and circles; the arrows indicate the generative relationship between these quantities.
By Bayes’ Theorem, the full posterior on both population parameters and object parameters is written
where is a prior on the population parameters; is the probability of being selected given the data of that object; is the likelihood of the data of an object given its object parameters; is the probability of object parameters given the population parameters; finally, is the normalization of , and depends on the selection function and the population parameters. When writing the data or object parameters without an index , we refer to the complete set of objects, .
4.1 Object parameters
Each WD has a set of object parameters encapsulated in , as listed in Table 1. Because the angular position errors can be neglected, the spatial position only varies in terms of the object’s distance. Conceptually, the most straightforward parametrization would be to have the distance as an object parameter. However, this creates some sampling difficulties arising from the fact that and are highly degenerate, especially when there is no parallax information available. In this work, we sample the object parameter posteriors using a Metropolis-Hastings Markov Chain Monte-Carlo algorithm (Metropolis et al., 1953; Brooks et al., 2011), which is more efficient when the posterior closer to a multivariate Gaussian in shape. This can be accomplished by a coordinate transformation, as is illustrated in Fig. 4. In the upper panel, the distance is parametrized in terms of , which is the relative change to the ideal distance given and . Let be the distance that maximizes the colour factor of the likelihood function, which is
The function is the inverse of , the difference between apparent and absolute magnitude. The distance parametrized by the object parameters is . It is crucial to account for the Jacobian factor that arises with this parametrization, in which the differentials are replaced according to
where the Jacobian is
which has no dependence on .
4.2 Object posterior
The posterior on a specific object also includes the term , normalized by the quantity . It is given by the population parametrized by
where is the number density of WDs as a function of spatial position. It is implicit in this expression that the true parallax is a function of the object parameters .
The normalization factor is given by an integral (or sum, in the case of a discrete variable) over the possible properties of a hypothetical WD drawn from the population model, multiplied by the probability of being selected, as
The selection function, , is the probability of being included in the sample, given by an object’s intrinsic properties and the sample construction cuts on observables.
5 Mock data and analysis
To test the algorithm, mock data are generated from the population model. While the exact values of the population parameters are of lesser importance (the main focus being the statistical method), we chose values with the following motivations.
For the temperature distribution we chose . To first order, WDs cool at a rate of , according to Mestel (1952), although numerical models differ from this cooling rate especially for cooler WDs.
For the distribution of surface gravity, we chose , , and . The surface gravity distribution mean and dispersion ( and ), motivated by Eisenstein et al. (2006). The tails of the surface gravity distribution are observed to be heavier than those of a Gaussian distribution, motivating a value .
For the fraction of DB WDs, we chose . The fraction of DB WDs is observed to vary with temperature, roughly in the range of 10–20 per cent (Bergeron et al., 2011).
We compare the case of having no astrometric information, versus having parallax measurements with the precision of Gaia DR2.
5.1 Sample selection
We define a sample by making cuts on observable properties, in our case on mock data of a generated catalogue. Depending on the errors of these observables, these cuts correspond more or less well to cuts in the objects’ intrinsic properties. We do not make a volume-limited sample by making cuts on parallax – we wish to compare with the case where astrometry is not available, hence we need to construct the sample without such information. We make cuts in observed apparent magnitude and observed colour. The cuts in colour correspond to upper and lower limits on the temperature of WDs included in our sample. The limit in apparent magnitude sets a maximum distance for a WD, as a function of temperature, surface gravity and type.
There are several reasons for not allowing very high temperatures in our samples (although the exact limit is rather arbitrary). Very hot objects are rare in terms of spatial density, but because they are seen to much greater distances they can still be numerous in a sample that is not volume-limited. How many are seen depends on the properties of the population, but this is degenerate with the assumed spatial distribution and the distribution of Galactic dust. Furthermore, with sufficiently high temperature, the peak of an object’s spectrum is of shorter wavelength than the band, in which case the inference on an object’s temperature becomes very weak. When working with actual data, it is also necessary to identify what objects are WDs and what objects are not. With very hot, far away objects this is difficult, especially since the distance will be poorly constrained. These issues can be circumvented with good parallax measurements, enabling the construction of a volume-limited sample.
We make the following cut in colour, demanding that
Without measurement errors, this cut corresponds to limiting the temperature of a DA type WD to K; for a DB WD the upper limit in temperature is less restrictive, as can be seen in Fig. 5.
We also make a cut in brightness, by demanding that the Gaia band apparent magnitude fulfils . In principle, this criterion could equally well have been formulated in terms of some combination of the apparent magnitudes.
Given these cuts on observables, the selection function is
where and are true observables with an implicit dependence on the object parameters. The error on is given by
assuming that the different magnitude errors are uncorrelated.
5.2 Generating a mock sample
The mock sample of WDs is generated by rejection sampling. An object is drawn randomly from the true population model: the object parameters , and are distributed as described in Eq. (2) and Eq. (5) and can be randomized analytically; the position is distributed according to and is randomized by rejection sampling, knowing that there is a maximal distance a WD can have in order to be included in our sample (observational errors included). A randomly drawn object is then assigned observable quantities, with errors as described in Sec. 3. If the object observables fulfil the selection cuts, the object is included in the sample; if not, it is rejected.
We construct a sample with 10,000 WDs. The distribution of true object parameter values is shown in Fig. 6, where selection effects are manifest. The maximum distance is clearly seen as a function of temperature, where hotter objects are seen further away. Due to the colour cut, the high temperature tail is more pronounced for the DB sub-population. It is also clear that low mass WDs are more likely to be included as they are more luminous, giving rise to some asymmetry in the distribution of surface gravity. The hottest object in this sample has an effective temperature of K and is of DB type. The most distant object is located at 1.56 kpc, is of DA type and has an effective temperature K surface gravity .
We infer the population parameters of a Bayesian hierarchical model, as described in Sec. 2, for our mock data sample, using a Monte-Carlo Markov Chain (MCMC) to trace the posterior distribution. We use a special purpose sampling algorithm called Metropolis-within-Gibbs (Gelman et al., 2013), in which the population parameters () and the object parameters of each individual WD () are varied separately, producing a computationally efficient marginalization of the object parameters. The WD type (), which is a discrete variable, alternates at random for each attempted step of the object parameter Metropolis step. For all other variables, which are continuous, the step is drawn at random from a multivariate normal distribution. The covariance matrices of these multivariate Gaussians, one for the population parameters and one for each of the 10,000 objects, are tuned in a thorough burn-in phase. For each iteration of the algorithm, the algorithm attempts 40 steps in population parameter space, and 40 steps in the parameter space of each separate object.
The population parameter prior, , is taken to be uniform and wide in all parameters around the true parameter values, with the exception of the lower bound to which is set to 1.
In both cases, with and without parallax information, the correct population parameter values are recovered by the posterior distribution. The inference on quantities (parametrizing the distribution of effective temperature) and (the fraction of DB WDs) is not significantly affected by also adding parallax information. For the remaining parameters, , , and (parametrizing the distribution of surface gravity), the inference is strongly affected when including parallax information. The width of the posterior distribution is roughly halved for all three quantities. There is a clear degeneracy between and , following the relation where the total variance on the surface gravity, , is preserved.
In Table 2, we present highest posterior density (HPD) intervals for the population parameters of three mock data samples, the first one being the sample discussed above. There is no indication of any bias in the inference of the population parameters.
6 White dwarf subpopulations
The population of WDs in the Milky Way is richer and more complicated than the model described in Sec. 2. While we do consider WD sub-populations in the sense of accounting for the difference between DA and DB type WDs, there are other meaningful ways to construct the population model. While we assume that the distribution of temperatures and surface gravities is the same between DA and DB WDs, derived from the same population parameters, it could be meaningful to have separate sets of population parameters for the two types.
In the same vein, one would expect the disk and halo WD population to have different properties. Furthermore, there is expected to be a sub-population of binary WD systems. We discuss how to model these two cases below.
6.1 Disk and halo populations
It would be interesting to see qualitative differences between disk and halo WDs. For example, the kinematically warmer halo population has older stars. The population of very old WDs is scientifically interesting, as it holds information about the star formation and dynamical history of the Milky Way.
In this population model, each sub-population would have its own set of population parameters: . They would each have their respective spatial number density distributions: and . It would be necessary to have a population parameter that describes the relative number density fraction of the two sub-populations at some reference point (for example the Sun’s position), such that they can be normalized. The posterior, analogous to Eq. (9) would read
The total number count is simply the sum over the two sub-populations, like
6.2 Binary population
Unresolved binary WD systems can be identified using only photometry and astrometry, in a similar way to the method presented in Widmark et al. (2018). For a binary system, the likelihood is the same as in Eq. (8), the difference being that the ugriz apparent magnitudes of the two component stars are added together, according to
where and are the -band apparent magnitudes of the two component stars.
The posterior density of a binary system will be written in terms of 7 parameters instead of 4, as we have temperature, surface gravity, and type of the two component stars, and the distance of the binary system.
We construct a population of mock binaries by random pairing of the singles population and the same selection criteria, although we also add a constraint in terms of cooling time of the two component stars. In addition to the effective temperature and surface gravity distributions of the two component stars, as described by Eq. (2), the probability of pairing also has a factor
where is the cooling time of the component WD (and equivalently for component ), as given by the Bergeron atmospheric model. The chosen time difference of million years prohibits the pairing of extremely cool and faint WDs with hotter ones. This is a reasonable assumption, as binary stars are typically born in the same system and with similar properties. Without this constraint, one component star would almost always be extremely faint, making binary identification almost impossible. For reference, the cooling time of a WD with and is roughly 500 Myr.
Figure 9 shows a receiver operating characteristic (ROC) curve for identification of binaries, for the cases of having and not having parallax information. The binaries are inferred with knowledge of the underlying population model, in the sense that the population parameters are known. The Bayesian evidence for being a single WD and being a binary WD is computed by sampling the posterior over the object parameters with an MCMC, and then approximating the integral by a first order multivariate Gaussian approximation, given by the covariance matrix and maximum posterior value MCMC chain. The integral is computed for all possible DA and DB combinations separately, with a multiplicity constraint for binary WDs, circumventing issues of multimodal posterior densities. This is done for 1000 mock data single WDs and 1000 mock data binary WD systems. It is clear from Fig. 9 that binary identification is significantly improved with parallax information, for which some binaries can be strongly identified even with a very low contamination rate (20 per cent binary identification with 0.1 per cent contamination).
This identification is made on mock data when the underlying population model is known. Working with real data brings many complications, not least coming from the fact that the population model is unknown and inferred. Even so, this test shows that it should be possible to identify a WD binary population. An important aspect that is not accounted for here is that some WD binaries are not drawn from the same distribution of masses as the population of single WDs. A tight binary system goes through phases of mass transfer and shared envelopes; thus there will be binaries with component WD with low mass and surface gravity. Such a binary system would actually be even easier to detect using this method, as they would be brighter due to multiplicity, and brighter still from being low mass and larger in size.
In this paper, we demonstrate how to infer properties of the local WD population using only astrometric and photometric information, in the framework of a Bayesian hierarchical model.
In our mock sample, we have limited ourselves to a total number of 10,000 WDs and a simple population model, in order to demonstrate the statistical method. The catalogue of WDs in a Gaia and SDSS cross-matched sample is expected to be around an order of magnitude larger, enabling us to fit a significantly more complicated model. The model could be extended with more complex distributions of effective temperature, surface gravity, and type, and by including sub-populations as discussed in Sec. 6. With a kinematic model, proper motion information would be very informative, especially in terms of differentiating between disk and halo WDs.
When working with real data, there are complications that are not included here but would be straightforward to implement within this framework. Most WD seen by Gaia and SDSS are very close to the Sun and almost unaffected by dust. However, hotter and more luminous WDs are seen to further distance and subject to dust reddening and extinction. With a good dust map, selection effects and photometric reddening for such objects can be accounted for. Also not included in this work are incompleteness effects, which are severe for WDs in Gaia DR2. This will improve significantly with future data releases, but will still be crucial to account for.
Gaia parallax measurements provide robust identification of WDs, enabling the construction of volume-limited samples, and breaks the degeneracy between distance and size. It is possible to differentiate sub-populations of WDs using this method, such as a population of binary WD systems. Our statistical model fully and correctly accounts for selection effects and observational uncertainties, permitting the construction of a large data sample, without the need to exclude objects with low signal-to-noise or missing parallax information.
We would like to thank Pierre Bergeron, for providing WD atmospheric models with SDSS and Gaia passband photometry. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. This work was also partially supported by a grant from the Simons Foundation.
- Bergeron et al. (1995) Bergeron P., Wesemael F., Beauchamp A., 1995, Publ. Astron. Soc. Pac., 107, 1047
- Bergeron et al. (2001) Bergeron P., Leggett S. K., Ruiz M. T., 2001, Astrophys. J. Suppl., 133, 413
- Bergeron et al. (2011) Bergeron P., et al., 2011, ApJ, 737, 28
- Brooks et al. (2011) Brooks S., Gelman A., Jones G., Meng X.-L., 2011, Handbook of Markov Chain Monte Carlo. CRC press
- Carrasco et al. (2014) Carrasco J. M., Catalán S., Jordi C., Tremblay P.-E., Napiwotzki R., Luri X., Robin A. C., Kowalski P. M., 2014, A&A, 565, A11
- Dame et al. (2016) Dame K., Gianninas A., Kilic M., Munn J. A., Brown W. R., Williams K. A., von Hippel T., Harris H. C., 2016, MNRAS, 463, 2453
- Eisenstein et al. (2006) Eisenstein D. J., et al., 2006, ApJS, 167, 40
- El-Badry et al. (2018) El-Badry K., Rix H.-W., Weisz D. R., 2018, preprint, (arXiv:1805.05849)
- Finley et al. (1997) Finley D. S., Koester D., Basri G., 1997, Astrophys. J., 488, 375
- Fontaine et al. (2001) Fontaine G., Brassard P., Bergeron P., 2001, PASP, 113, 409
- García-Berro & Oswalt (2016) García-Berro E., Oswalt T. D., 2016, New Astron. Rev., 72, 1
- Gelman et al. (2013) Gelman A., Carlin J., Stern H., Dunson D., Vehtari A., Rubin D., 2013, Bayesian Data Analysis: Third Edition
- Gentile Fusillo et al. (2018) Gentile Fusillo N. P., et al., 2018, preprint, (arXiv:1807.03315)
- Isern et al. (1998) Isern J., García-Berro E., Hernanz M., Mochkovitch R., Torres S., 1998, ApJ, 503, 239
- Jordan (2007) Jordan S., 2007, ASP Conf. Ser., 372, 139
- Jurić et al. (2008) Jurić M., et al., 2008, ApJ, 673, 864
- Kepler et al. (2015) Kepler S. O., et al., 2015, MNRAS, 446, 4078
- Kleinman et al. (2013) Kleinman S. J., et al., 2013, ApJS, 204, 5
- Koester & Kepler (2015) Koester D., Kepler S. O., 2015, A&A, 583, A86
- Kruijssen et al. (2018) Kruijssen J. M. D., Pfeffer J. L., Reina-Campos M., Crain R. A., Bastian N., 2018, MNRAS,
- Lindegren et al. (2018) Lindegren L., et al., 2018, preprint, (arXiv:1804.09366)
- Livio & Mazzali (2018) Livio M., Mazzali P., 2018, Phys. Rept., 736, 1
- Mestel (1952) Mestel L., 1952, MNRAS, 112, 583
- Metropolis et al. (1953) Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., 1953, J. Chem. Phys., 21, 1087
- Mortlock et al. (2009) Mortlock D. J., Peiris H. V., Ivezic Z., 2009, Mon. Not. Roy. Astron. Soc., 399, 699
- Perryman et al. (1997) Perryman M. A. C., et al., 1997, A&A, 323, L49
- Postnov & Yungelson (2014) Postnov K. A., Yungelson L. R., 2014, Living Rev. Rel., 17, 3
- Ritossa et al. (1996) Ritossa C., Garcia-Berro E., Iben Jr. I., 1996, ApJ, 460, 489
- Smartt et al. (2009) Smartt S. J., Eldridge J. J., Crockett R. M., Maund J. R., 2009, MNRAS, 395, 1409
- Tremblay & Bergeron (2008) Tremblay P. E., Bergeron P., 2008, Astrophys. J., 672, 1144
- Widmark et al. (2018) Widmark A., Leistedt B., Hogg D. W., 2018, ApJ, 857, 114
- Winget et al. (1987) Winget D. E., Hansen C. J., Liebert J., van Horn H. M., Fontaine G., Nather R. E., Kepler S. O., Lamb D. Q., 1987, ApJ, 315, L77