<figure><img src="image_1.png" /></figure>

HAL
open science

# The Canada-France Imaging Survey: Reconstructing the Milky Way Star Formation History from Its White Dwarf Population

Nicholas Fantin, Patrick Côté, Alan Mcconnachie, Pierre Bergeron, Jean-Charles Cuillandre, Stephen Gwyn, Rodrigo Ibata, Guillaume Thomas, Raymond Carlberg, Sébastien Fabbro, et al.

► To cite this version:

Nicholas Fantin, Patrick Côté, Alan Mcconnachie, Pierre Bergeron, Jean-Charles Cuillandre, et al. The Canada-France Imaging Survey: Reconstructing the Milky Way Star Formation History from Its White Dwarf Population. The Astrophysical Journal, 2019, 887 (2), pp.148. 10.3847/1538-4357/ab5521. hal-03013810

HAL Id: hal-03013810

https://hal.science/hal-03013810v1

Submitted on 11 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The Canada-France Imaging Survey: Reconstructing the Milky Way Star Formation History from Its White Dwarf Population

Nicholas J. Fantin$^{1,2}$, Patrick Côté$^2$, Alan W. McConnachie$^2$, Pierre Bergeron$^3$, Jean-Charles Guillardre$^4$,

Stephen D. J. Gwyn$^2$, Rodrigo A. Ibata$^5$, Guillaume F. Thomas$^2$, Raymond G. Carlberg$^6$, Sébastien Fabbro$^2$, Misha Haywood$^7$,

Ariane Lançon$^8$, Geraint F. Lewis$^9$, Khyati Malhan$^{10}$, Nicolas F. Martin$^{11,12}$, Julio F. Navarro$^1$, Douglas Scott$^{13}$, and

Else Starkenburg$^{14}$

$^1$Department of Physics and Astronomy, University of Victoria, Victoria, BC V8P 1A1, Canada; nfantin@uvic.ca

$^2$National Research Council of Canada, Herzberg Astronomy & Astrophysics Research Centre, 5071 W. Saanich Road, Victoria, BC V9E 2E7, Canada

$^3$Département de Physique, Université de Montréal, C.P. 6128, Succ. Centre-Ville, Montréal, QC H3C 3J7, Canada

$^4$AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, Observatoire de Paris, PSL University, F-91191 Gif-sur-Yvette Cedex, France

$^5$Observatoire de Strasbourg, CNRS, UMR 7550, F-67000 Strasbourg, France

$^6$Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada

$^7$GEPI, Observatoire de Paris, PSL Research University, CNRS, Place Jules Janssen, F-92190 Meudon, France

$^8$Observatoire Astronomique, 11 rue de l'Université, F-67000 Strasbourg, France

$^9$Sydney Institute for Astronomy, School of Physics, A28, The University of Sydney, Sydney, NSW 2006, Australia

$^{10}$The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden

$^{11}$Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France

$^{12}$Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

$^{13}$Department of Physics and Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada

$^{14}$Leibniz Institute for Astrophysics Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany

Received 2019 August 19; revised 2019 November 4; accepted 2019 November 5; published 2019 December 17

## Abstract

As the remnants of stars with initial masses $\lesssim 8 M_\odot$, white dwarfs contain valuable information on the formation histories of stellar populations. In this paper, we use deep, high-quality, u-band photometry from the Canada-France Imaging Survey, griz photometry from Pan-STARRS1, as well as proper motions from Gaia DR2, to select 25,156 white dwarf candidates over $\sim 4500$ deg$^2$ using a reduced proper motion diagram. We develop a new white dwarf population synthesis code that returns mock observations of the Galactic field white dwarf population for a given star formation history, while simultaneously taking into account the geometry of the Milky Way (MW), survey parameters, and selection effects. We use this model to derive the star formation histories of the thin disk, thick disk, and stellar halo. Our results show that the MW disk began forming stars ($11.3 \pm 0.5$) Gyr ago, with a peak rate of ($8.8 \pm 1.4$) $M_\odot$ yr$^{-1}$ at ($9.8 \pm 0.4$) Gyr, before a slow decline to a constant rate until the present day—consistent with recent results suggesting a merging event with a satellite galaxy. Studying the residuals between the data and best-fit model shows evidence for a slight increase in star formation over the past 3 Gyr. We fit the local fraction of helium-atmosphere white dwarfs to be ($21 \pm 3$)%. Incorporating this methodology with data from future wide-field surveys such as the Large Synoptic Survey Telescope, Euclid, The Cosmological Advanced Survey Telescope for Optical and ultraviolet Research, and the Wide Field Infrared Survey Telescope should provide an unprecedented view into the formation of the MW at its earliest epoch through its white dwarfs.

Unified Astronomy Thesaurus concepts: White dwarf stars (1799); Milky Way Galaxy (1059); Milky Way evolution (1052)

## 1. Introduction

The formation and evolution of disk galaxies, and the Milky Way (MW) in particular, has long been an important topic in astronomy because these galaxies dominate star formation activity in the low-redshift universe. Previous investigations into the formation and evolution of MW-like galaxies have focused on the assembled mass as a function of redshift (see, e.g., van Dokkum et al. 2013), or compared resolved stellar populations to theoretical isochrones in local galaxies such as Andromeda (Ferguson et al. 2005; Brown et al. 2006). Within the MW itself, much of the attention has been on the stellar metallicity distribution because this, once combined with a model for the gas infall rate, contains information on the star formation rate (SFR) as a function of time (see, e.g., Snaith et al. 2015; Haywood et al. 2018; Toyouchi & Chiba 2018). Results from such Galactic chemical evolution (GCE) models typically show a strong period of initial star formation followed

by a slowly decaying trend toward the present day, features usually associated with a thick and thin disk, respectively.

The discovery of the dichotomy of the disk began with Gilmore & Reid (1983) who used star counts toward the South Galactic Pole and found that the vertical distribution was well described by two separate exponential profiles with different scale heights. Follow-up studies have shown that the thick disk stars are more alpha-rich, iron-poor, and have a larger velocity dispersion—evidence that they formed much earlier than their thin disk counterparts (Norris et al. 1985; Eggen 1998). Within the past decade, a renewed emphasis has been placed on the early formation of the MW’s disk. The disk dichotomy has been questioned by Bovy et al. (2012) who used mono-abundance populations to show that each population could be described by a single exponential scale height. Haywood et al. (2016), however, argued that the dip in the $[\alpha/Fe]$ distribution seen in the Apache Point Observatory Galaxy Evolution Experiment (APOGEE) can only be produced by two separate
epochs of star formation, and associates each epoch as belonging to the thin and thick disks, respectively (Snaith et al. 2015; Haywood et al. 2018). Haywood et al. (2016) further argues that these two results are consistent if one assumes that the scale height of the thick disk decreases as a function of time so that the differences in structural parameters between the young thick disk and old thin disk are minor. Such a decrease in the scale height, from 800 to 340 pc, was found to match star counts in the Sloan Digital Sky Survey (SDSS) and Two Micron All Sky Survey by Robin et al. (2014), with a mean scale height of 535 ± 50 pc. Further emphasis has been placed on this transitional epoch with recent results by Helmi et al. (2018) and Belokurov et al. (2018), who show that this epoch of the MW may have played host to a merger with at least one satellite galaxy, a scenario that would likely play an important role in the formation of the thick disk.

With a complete star formation history (SFH), the stellar mass of each component can be estimated. Traditionally, star count models calibrated to the solar neighborhood have found that the thick disk contributes about 12% of the local mass (Jurić et al. 2008). However, a recent study by Snaith et al. (2015) using the chemical abundances of local F, G, and K-type stars concluded that the thin and thick disk masses may be comparable. This discrepancy likely arises due to the regions of the MW that are probed by each survey. Many models use data confined to the solar neighborhood, a location at which the thin disk dominates, and thus extrapolating local results to the entire disk becomes problematic.

While the SFH is imprinted in the chemical enrichment of the MW, stars also leave behind another important tool for astronomers, white dwarfs. White dwarfs form at the end of the asymptotic giant phase as the outer layers are ejected to expose the degenerate carbon-oxygen core. Because these stars are remnants of stars with initial masses less than 8 $M_\odot$, they contain valuable information regarding the formation history of stellar populations as a whole.

Because their evolution is well understood—the age is primarily a function of temperature, luminosity, and atmospheric composition—white dwarfs have been used in the past to estimate ages. These techniques have been applied to various stellar populations, from globular clusters (see, e.g., Richer et al. 2006; Hansen et al. 2007; Bedin et al. 2009) to MW components (Kilic et al. 2017). The first such study of the MW’s disk was carried out by Winget et al. (1987) who used the Luyten Half-second Proper Motion Catalog to estimate an age of (9.3 ± 2.0) Gyr for the Galactic disk, and (10.3 ± 2.2) Gyr for the age of the universe, where the age refers to the onset of star formation. These early studies, however, suffered from strong incompleteness and small number statistics due to shallow photometric data. A more recent analysis by Kilic et al. (2017) used the white dwarf luminosity function (WDLF) from Munn et al. (2017) to simultaneously derive the ages of the MW components, obtaining age measurements of (7.8 ± 0.4) Gyr for the thin disk and (9.7 ± 0.2) Gyr for the thick disk. Due to the small number of halo white dwarfs in their sample, the age of the halo was essentially unconstrained.

With a large enough sample, the complete SFH of a stellar population can be determined from white dwarfs in two ways. First, with accurate distances one can calculate the luminosity, radius, and effective temperature and convert this into a cooling age (Giammichele et al. 2012). The total age of a white dwarf

can be found by adding the cooling age to the progenitor age, which is found by converting the white dwarf mass to a progenitor mass using the initial-to-final mass relation (IFMR). Tremblay et al. (2014) applied this method to the local (20 pc) sample and measured the SFH in the solar neighborhood. Their results suggest a roughly constant SFR over the past 10 Gyr, with a slight drop 5 Gyr ago (see their Figure 9). This method, however, is typically constrained to small volumes given the distance uncertainties.

The second method is to invert the WDLF. Noh & Scalo (1990) showed that varying the functional form for the SFR would leave imprints on the resulting WDLF, such as the slope and turn-off. This method was used by Rowell (2013) to invert the SDSS WDLF assembled by Harris et al. (2006), and showed a resulting SFH composed of two broad peaks (at 2 and 9 Gyr) separated by a significant dip. The advantage of using this method is that the volume that can be probed is considerably larger; the downside, however, is that volume and completeness corrections must be applied. The recent advances made in wide-field astronomy present an opportunity to increase the sampled volume, and the combination of deep, multiband photometry and accurate proper motions are ideal for selecting white dwarfs.

In this paper, we introduce a third technique: forward modeling. This involves simulating a sample of white dwarfs given a star formation prescription. The advantage of forward modeling is that complex systematics, like completeness, observational uncertainties, selection effects, and the structure of the MW can all be accounted for. We use data from the Canada-France Imaging Survey (CFIS), Pan-STARRS1, and Gaia DR2 to derive the SFH of the MW components using a new population synthesis code. Section 2 describes the data set, while Section 3 describes our white dwarf population synthesis code. We present the results in Section 4, discuss their implications in Section 5, and summarize in Section 6.

## 2. Data

The photometric data used in this paper were acquired as part of the CFIS (Ibata et al. 2017) and the *grizy* Pan-STARRS1 (PS1) 3$\pi$ survey (Chambers et al. 2016). CFIS is an ongoing large program at the Canada-France-Hawaii Telescope (CFHT) that aims to obtain 10,000 deg$^2$ of u- and 5000 deg$^2$ of r-band photometry to $5\sigma$ point source depths of 24.2 and 24.85 AB mag, respectively. The survey is performed under excellent conditions, with a median u-band image quality of 0.78. This paper makes use of the CFIS-u catalog, which covers ~4500 deg$^2$ as of the end of the 2018A semester. The footprint is highlighted in gray, in both equatorial and Galactic coordinates, as part of Figure 1.

The CFIS-u catalog was merged with PS1 *grizy* forced photometry as described in Thomas et al. (2018). Given that the PS1 and *Gaia* catalogs cover the whole sky as seen from the CFHT, our area is set by the CFIS-u footprint. This catalog retains more than 98% of PS1 detections within the magnitude range of interest. Stars are selected using the recommended star-galaxy classification from Farrow et al. (2014): $i_{PSF} - i_{Kron} < 0.05$ mag.

We also merge our catalog with *Gaia* DR2 (Gaia Collaboration et al. 2018) in order to obtain proper motions. We apply the same astrometric quality cuts on the data as those
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left: equatorial positions of thin disk stars within our model. Right: the model in Galactic coordinates. In both panels, the shaded region represents the CFIS-*u* footprint as of the end of the 2018A semester.</figcaption></figure>

presented in Gentile Fusillo et al. (2019),

$$ \begin{array}{l} \text{ASTROMETRIC\_SIGMA5D\_MAX} < 1.5 \text{ OR} \\ \text{(ASTROMETRIC\_EXCESS\_NOISE} < 1 \\ \text{AND PARALLAX\_OVER\_ERROR} > 4 \\ \text{AND SQRT(PMRA}^2 + \text{PMDEC}^2) > 10 \text{ mas,} \end{array} \qquad (1) $$

which removes objects with poor astrometric measurements. Specifically, ASTROMETRIC\_SIGMA5D\_MAX represents the longest axis of the 5d error ellipsoid, in mas, and large values indicate that one of the astrometric parameters is poorly constrained. ASTROMETRIC\_EXCESS\_NOISE is a measure of the difference between the best-fit astrometric solution and the data, with large values suggesting a poor astrometric solution. Finally, PARALLAX\_OVER\_ERROR is the relation between the parallax measurement and its associated error. The limiting magnitude of Gaia (G = 20.7) represents the primary limitation of our catalog, and the implications of this are discussed in Section 5. Our CFIS-PS1-Gaia catalog contains more than 21.5 million unique point sources down to the Gaia limiting magnitude.

## 2.1. White Dwarf Selection

Using our combined catalog, we have constructed a reduced proper motion diagram (RPMD) based on photometry from CFIS and PS1, and proper motions from Gaia DR2, to separate white dwarfs from other point sources. The reduced proper motion, H (Luyten 1922), combines the apparent magnitude, m, and the proper motion, µ, and can be used as a proxy for absolute magnitude, M, and tangential velocity, vᵢ.

$$ \begin{array}{l} H = m + 5 \log \mu + 5 \\ = M + 5 \log v_t - 3.379. \end{array} \qquad (2) $$

Due to their intrinsic faintness, for a given temperature white dwarfs are well separated in the RPMD from main-sequence stars with similar colors (see, e.g., Harris et al. 2006; Rowell & Hambly 2011; Fantin et al. 2017; Munn et al. 2017). Using model cooling tracks we select white dwarfs with tangential velocities greater than 20 km s⁻¹ as shown in the left-hand panel of Figure 2. These synthetic magnitudes have been calculated in the CFIS-*u* and PS1 grizy bands for H- and He-atmosphere models using the procedure described in Holberg & Bergeron (2006). The white dwarf cooling sequences are similar to those described in Fontaine et al. (2001) with (50/50) C/O-core compositions, M<sub>He</sub>/M<sub>⋆</sub> = 10⁻², and M<sub>H</sub>/M<sub>⋆</sub> = 10⁻⁴ or 10⁻¹⁰ for H- and He-atmosphere

white dwarfs, respectively.<fnref n="15" /> Our selection procedure identifies 25,156 white dwarf candidates within the CFIS-*u* footprint. The resulting sample, along with model H- and He-atmosphere tracks, can be seen in color–color space in the right-hand panel of Figure 2.

Previous selections using the RPMD have obtained contamination fractions from non-white dwarfs (typically hot sdO, sdB, or cooler main-sequence stars) of a few percent (e.g., Harris et al. 2006) at a velocity cut of vᵢ > 30 km s⁻¹, with the rate increasing for lower velocity cuts (Kilic et al. 2010). Given that Gaia DR2 has much better astrometry than previous studies, we have chosen to lower our velocity cut to vᵢ > 20 km s⁻¹. Matching our objects to SIMBAD, we find 3688 spectroscopic matches, of which 3589 are white dwarfs. The remaining objects are a mix of hot subdwarfs (69 objects), QSOs with spurious proper motion measurements (20 objects), and 10 miscellaneous objects, including main-sequence stars. The resulting contamination rate is thus at least ~3%, a value that has been incorporated into the results presented in Section 4.

The color–color diagram highlights the precise photometry of CFIS and PS1 and the ability to separate DA (H-lines present) and DB (He-lines present) white dwarfs between −0.6 < g<sub>PS</sub> − i<sub>PS</sub> < 0.0 mag. The RPMD also highlights the precise astrometry of Gaia DR2, which results in a clean sequence of white dwarfs in color–color space. In Figure 3 we present the observational uncertainties in both the photometry and astrometry as a function of magnitude for our white dwarf sample, highlighting the depth of the CFIS-*u* band and the precision of the Gaia DR2 proper motions.

## 2.2. Completeness

We calculate the completeness within our data set by comparing our CFIS-PS1-Gaia catalog to the CFIS-PS1 catalog. Thomas et al. (2018) showed the PS1 g-band is complete to 22.0 mag, which is deeper than our sample (see Figure 3) and thus we assume that the CFIS-PS1 data set is complete over the magnitude range covered by the Gaia survey. Our completeness curves for the PS1 g and i bands, as well as the CFIS-*u* band, can be seen in Figure 4.

## 3. Modeling the MW’s White Dwarf Population

While the MW’s white dwarf population contains important information regarding its evolutionary history, it is generally

<footnote n="15">¹⁵ See http://www.astro.umontreal.ca/~bergeron/CoolingModels.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: CFIS-PS1-Gaia RPMD showing all 24.5 million sources. White dwarfs for this study were selected if they had inferred tangential velocities greater than 20 km s⁻¹ (the region bound by the blue line). Right: color–color diagram showing the resulting white dwarf candidates selected from the RPMD. For reference, 0.6 M☉ H- and He-atmosphere model tracks are plotted in red and cyan, respectively.</figcaption></figure>

<figure><img src="image_6.png" /></figure>

not the focus of stellar population synthesis codes. Prominent population synthesis codes, like the Besançon model (Robin et al. 2003, 2014), add white dwarfs given an observed local density, while TRILEGAL model (Girardi et al. 2005) has issues reproducing observations of white dwarfs (Fantin et al. 2017). Thus, a model that allows white dwarfs to form via the evolution of progenitor stars, as opposed to being added manually, allows one to study various input parameters, such as the SFH, which can affect the resulting white dwarf population. In this section, we describe our white dwarf population synthesis code and detail our various input assumptions about the MW.

### 3.1. Model Functions

We construct a synthetic model of the MW’s white dwarf population in order to compare to our data set. We start by assuming a three-component Galaxy consisting of a thin disk, a thick disk, and a stellar halo, with density distributions given in Table 1. We assume a double exponential profile in both R and z for the thin and thick disk. We assume a scale length of 2.3 kpc for both the thin and thick disk, and scale heights of 300 and 550 pc, respectively, similar to the mean values presented by Robin et al. (2014). For each component, we generate stars at positions and masses given by the assumed density distributions and Kroupa (2001) initial-mass function (IMF).

Stars are spawned at Galactocentric radii (R) and heights above the plane (z) appropriate for the distributions shown in Table 1. The stars are evenly distributed in the angular coordinate, with 0° representing the line connecting the solar position in the plane to the Galactic center. The Galactocentric coordinates are then converted to R.A. and decl., and a distance from the Sun is calculated assuming the Sun is at a position of (R, θ, z) = (8340 pc, 0 pc, 17 pc) (Reid et al. 2014; Karim & Mamajek 2017). Figure 1 shows the resulting equatorial (left)

and Galactic (right) positions for the thin disk with the CFIS-u footprint highlighted in gray.

The amount of time it will take each star to become a white dwarf, the progenitor lifetime, is calculated using the analytic stellar lifetimes from Hurley et al. (2000). This functional form takes the initial mass and the metallicity of each star and returns the lifetime of the star. For our model, we assume solar metallicity for the thin disk, [Fe/H] = −0.7 for the thick disk, and [Fe/H] = −1.5 for the halo (Peng et al. 2013). Each star is also assigned a birth date (its formation age) that is randomly generated given the functional form for the SFR. We assume a skewed Gaussian SFH for each component,

$$SFR(t) = \rho_0 \frac{2}{\sigma_t} \phi\left(\frac{t-\xi}{\sigma_t}\right) \Phi\left(\alpha\left(\frac{t-\xi}{\sigma_t}\right)\right), \quad (3)$$

where φ(t) is the standard normal probability distribution that is symmetric about the mean, ξ, and Φ(t) is its cumulative distribution function. The skewness parameter is α, the standard deviation is given by σt, and ρ₀ is the space density.

A star will have formed a white dwarf if the lifetime up to the end of the AGB phase is less than its formation age, with this difference being the white dwarf cooling time. The white dwarfs are then given space velocities according to the distributions presented in Table 1. These solar position values were adopted from Robin et al. (2017), who derived them by comparing the Besançon model to kinematic data from RAVE and Gaia DR1. These values are similar to those derived by Rowell & Hambly (2011) using white dwarfs in the Super-COSMOS survey. We do not allow the velocities to vary in R or z, as shown in Robin et al. (2017), because our sample is local (see Section 4) and therefore this gradient would have a minimal impact on our results. We also do not allow our thin disk velocity dispersion to vary with age, but instead take a mean value because the small change would have little effect on a star’s position in the RPMD. These velocities are
<figure><img src="image_7.png" /><figcaption>Figure 3. Observational uncertainties in (a) CFIS u, (b) PS1 g, (c) PS1 i, and (d) proper motions as a function of magnitude in our white dwarf sample. The model, as described in Section 3, samples a Gaussian at each point with a mean (given by the red line showing a polynomial fit) and standard deviation in order to obtain an uncertainty value for each mock white dwarf.</figcaption></figure>

converted to a proper motion using the equations from Johnson & Soderblom (1987).

The white dwarf mass, in solar masses, is calculated from the initial progenitor mass using the IFMR from Kalirai et al. (2008).

$$M_{WD} = (0.109 \pm 0.007)M_i + (0.394 \pm 0.025)M_{\odot}. \quad (4)$$

The final parameter needed to calculate the photometry of a white dwarf is its spectral type. We include pure H- and He-atmosphere white dwarfs in our model, and each white dwarf is designated as one or the other given an input fraction, $f_{He}$. With the cooling age, white dwarf mass, and type in hand, we determine the absolute magnitudes in the required bands using the white dwarf cooling models as described in Section 2 and seen in Figure 2. We then determine the reddening in each band at the white dwarf's position using the extinction coefficients and E(B-V) values from the Bayestar 3D dust maps of Green et al. (2015, 2018), before converting the model absolute

magnitudes to apparent magnitudes. We then add experimental uncertainties based on the observed relation between magnitude and error in the given photometric bands and proper motion (see Figure 3). Finally, completeness (see Figure 4) and selection effects (see Figure 2) are applied based on the observations.

The model returns a catalog of simulated white dwarfs within a given footprint. Catalog information includes positions, distances, masses, cooling ages, spectral types, proper motions, and apparent magnitudes in the given photometric bands. Typical results can be seen in Figure 5, which shows color-color diagrams (left column) and RPMDs (right column) for young (top row), intermediate-age (middle row), and old populations (bottom row).

While the above description has been tailored for this particular study, we note that the model can readily be adapted for any future studies by modifying the bandpass and survey parameters.
<figure><img src="image_8.png" /><figcaption>Figure 4. The resulting completeness as a function of magnitude for our CFIS-PS1-Gaia sample in the CFIS u (blue), PS1 g (green), and PS1 i bands (red), respectively. These completeness functions were calculated assuming that the CFIS-PS1 catalog is complete over the magnitude range of Gaia DR2, which was shown to be true in Thomas et al. (2018).</figcaption></figure>

## 3.2. Fitting Method: Approximate Bayesian Computation Markov Chain Monte Carlo (MCMC)

We begin by parameterizing the SFH for each Galactic component following Equation (3). This results in four parameters for each component: the mean, $\xi$, which we will call the functional age, the scale, $\sigma_t$, the skewness $\alpha$, and the SFR. This results in a total of 12 parameters. We also fit for the ratio between the hydrogen and helium-atmosphere white dwarfs, which we call the He fraction, $f_{\text{He}}$. This fraction is assumed to be equal for all three populations. This brings the total number of parameters to 13.

We integrate a given SFH to determine the total stellar mass and select the corresponding number of stars contained in each MW component. For each star, we then generate a formation age given the SFH. We then remove every star that does not fall within the CFIS footprint and/or will not have formed a white dwarf at the present day. The cooling age for those stars that remain is then calculated as being the difference between the formation age and the progenitor lifetime.

With the cooling ages, white dwarf masses, and spectral type, we then calculate apparent magnitudes in the CFIS and PS1 bands using the model cooling curves described in Section 2. Observational uncertainties are added given the relations shown in Figure 3. A photometric completeness correction is then made given the curves presented in Figure 4.

Finally, we perform the same selection methods on the model as we do on the data. We begin by applying identical distributions in the Gaia parameters included in Equation (1) to our model white dwarfs before making the same selection as presented in Equation (1). Using the proper motion we calculate the reduced proper motion using Equation (2) and select those stars whose tangential velocity is greater than the 20 km s⁻¹ curve (see Figure 2).

We use astroABC (Jennings & Madigan 2017), a python-based approximate Bayesian computation (ABC) MCMC approach in order to calibrate our model. This so-called “likelihood-free” method is useful in cases where a likelihood function is difficult to calculate, or unknown altogether. Unlike traditional MCMC methods used in astronomy (see, e.g., emcee; Foreman-Mackey et al. 2013), ABC MCMC does not calculate an explicit likelihood but instead relies on forward modeling. Given a set of model parameters, a simulated data set is produced and compared to the data in order to either accept or reject the proposed parameters. More specifically, if we draw parameters $\theta$, we accept these parameters if $\rho(D - D^*(\theta)) < \epsilon$, where $\rho$ is the distance metric between the data set, $D$, and the simulated data given $\theta$, $D^*(\theta)$, and $\epsilon$ is called the tolerance threshold. The distance metric used to compare the data and the model is

$$ \rho(D_i - D_i^*(\theta)) = \sum_i \frac{(D_i - D_i^*(\theta))^2}{2\sigma_i^2}. \quad (5) $$

The result at each step is a set of independent samples below a generated tolerance threshold (Marjoram et al. 2003).

The algorithm starts off with a large threshold and shrinks it at each step in order to approximate the PDF. When this threshold is exactly zero, then one is sampling directly from the posterior PDF; however, in that case, the acceptance fraction declines sharply as it is difficult to perfectly simulate a data set. On the other hand, if $\epsilon$ is too large, then the algorithm is sampling from the prior. Thus, the choice of $\epsilon$ determines how much of an approximation the resulting PDFs are while minimizing the required total computational time (Marjoram et al. 2003).

As described in Marjoram et al. (2003), ABC methods typically require lower-dimensional summary statistics, $S(D)$, which are representative of the whole data set in order to reduce computational costs. This summary statistic must encompass any required information contained within the data itself without neglecting anything important. Similar to Mor et al. (2018), we use a combination of binned colors, magnitudes, and proper motions. Specifically, we compare our real and simulated data using a 3D histogram composed of a color-color diagram ($g-i$, $u-g$), and the reduced proper motion, $H$. This combination of colors and proper motions was chosen in order to encompass both age and SFR information (colors), as well as the contributions from each component (apparent magnitudes and proper motions via $H$). Furthermore, the color-color diagram contains information regarding the He fraction in our sample, because they are well separated due to the exquisite $u$-band photometry provided by CFIS (see Figure 2). Thus, in Equation (5), $i$ represents each bin in the histogram, and $\sigma$ is the sum in quadrature of the uncertainty in both the model and data in each bin. For instance, where the model predicts a star, but we do not observe one, we use substitute $D_i+1$ into Equation (5) in order to penalize the model by increasing the distance. We then run astroABC with 1000 particles to obtain the approximate PDFs for the age, duration, and SFR for the three Galactic components, as well as the He fraction. We use uniform priors for every variable, with the only constraint being that the age of the components cannot exceed the age of the universe, which we set to 13.8 Gyr (Planck Collaboration et al. 2018), and that the scale parameters, $\sigma_t$, must be greater than zero.
<table><caption>**Table 1**

Assumed Model Distributions for Each Component</caption><thead><tr><th>Component</th><th>ρ(R, z)</th><th>ξ(M)</th><th>SFH</th><th>[Fe/H]</th><th>⟨V<sub>☉</sub>⟩<br>(km s<sup>-1</sup>)</th><th>σ<sub>U</sub><br>(km s<sup>-1</sup>)</th><th>σ<sub>V</sub><br>(km s<sup>-1</sup>)</th><th>σ<sub>W</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Thin disk</td><td>e<sup>-R/h<sub>R</sub></sup>e<sup>-z/h<sub>z</sub></sup><br>h<sub>R</sub> = 2300 pc, h<sub>z</sub> = 300 pc</td><td>Kroupa</td><td>Skewed Gaussian</td><td>0.0</td><td>-12</td><td>33</td><td>15</td><td>15</td></tr><tr><td>Thick disk</td><td>e<sup>-R/h<sub>R</sub></sup>e<sup>-z/h<sub>z</sub></sup><br>h<sub>R</sub> = 2300 pc, h<sub>z</sub> = 550 pc</td><td>Kroupa</td><td>Skewed Gaussian</td><td>-0.7</td><td>-85</td><td>40</td><td>32</td><td>28</td></tr><tr><td>Stellar halo</td><td>r<sup>-2.44</sup></td><td>Kroupa</td><td>Skewed Gaussian</td><td>-1.5</td><td>-226</td><td>131</td><td>106</td><td>85</td></tr></tbody></table>

Note. The solar position is assumed to be (R<sub>☉</sub>, θ<sub>☉</sub>, z<sub>☉</sub>) = (8340 pc, 0, 17 pc).

## 4. Results and Comparison to the Literature

The resulting corner plot containing the PDFs can be seen in Figure 6. The results from our fit can be seen in the bottom panels of Figure 7, which show a model instance with the mean parameters as estimated by the final step of the Markov chain (the central line in the histograms from Figure 6). The top two panels show the actual data for comparison. We are able to reproduce the thin and thick disk population quite well, although the halo parameters are much less constrained. This is a result of the *Gaia* limiting magnitude because the intrinsic faintness and scarcity of halo white dwarfs results in very few appearing in the data set. For reference, only 44 objects lie below the 200 km s<sup>-1</sup> model track in Figure 2, which is typically associated with the Galactic halo. This sample is smaller than that obtained with deeper SDSS photometry by (Munn et al. 2017), which was used by Kilic et al. (2017) to obtain a halo age of 12.5<sup>+1.4</sup><sub>-3.4</sub> Gyr. A larger sample will be required to better constrain the halo parameters—in particular, a sample that includes objects beyond the turnover in the halo luminosity function. This will require deeper proper motion surveys, such as PS1 DR2, which will be investigated in a future paper. Within the next decade, upcoming imaging surveys such as the *Wide Field Infrared Survey Telescope* (WFIRST; Spergel et al. 2015) and the Large Synoptic Survey Telescope (LSST; LSST Science Collaboration et al. 2009) will provide proper motions to unprecedented depths, allowing for a more comprehensive study of the Galactic halo.

### 4.1. *Our Sample, in the Context of the MW*

In order to properly compare with previous results, we begin by showing the volume that is sampled by our model. The resulting sight lines in both *R* (the distance from the center of the MW) and *z* (the distance from the midplane) can be seen in the left-hand panel of Figure 8. The cyan-dashed lines show heliocentric distances of 1, 2, and 3 kpc, respectively, showing that the vast majority of our sample lies within 1 kpc of the Sun. This is confirmed in the right-hand panel, which shows the distance distributions of our model white dwarfs. The result is that only 4% of white dwarfs in our sample lie at distances greater than 1 kpc, with a median distance of 388 pc. Given these values, we consider our sample to be representative of the solar neighborhood, and not necessarily the whole disk.

### 4.2. SFH

Our resulting SFH, integrated over the entirety of each component, can be seen in the left-hand panel of Figure 9, where we have sampled the resulting PDFs and plotted the mean (black) and 1σ contours (gray). Our results show an

initial period of thick disk star formation, lasting 3 Gyr and peaking at (9.8 ± 0.3) Gyr. This is followed by a nearly 1 Gyr decline in star formation before finishing with a nearly constant SFR for the thin disk at approximately 3.5 M<sub>☉</sub> yr<sup>-1</sup>. Below, we present a comparison between our SFH and those found in the literature. We make note of a caveat to these comparisons, which is that we assume that each component is defined by a unimodal SFH. Thus, if a given component contains multiple bursts, we will fit the average through the bursts (see Section 5 for further discussion).

The star formation, and in turn, the mass assembly, history of the MW has been an open question for decades. With the increase in accurate distances in the solar neighborhood from surveys like *Hipparcos* (ESA 1997), the local SFH could be derived by comparing the synthetic color–magnitude diagrams (CMD) to observations. This process was performed by Vergely et al. (2002) using *Hipparcos* stars brighter than V = 8, as well as Cignoni et al. (2006) who used all *Hipparcos* stars within 80 pc. Both of their samples are dominated by thin disk objects given their selection criteria, and in turn, they show a slowly rising star formation until it hits a peak 2 Gyr ago. These SFHs can be seen in the left-hand panel of Figure 10 as the red-dashed line and cyan-dotted line, respectively. In this figure, we have scaled our SFH and hence only the shapes should be compared. We note that neither SFH recovers a strong peak early-on in the formation of the MW (the Cignoni et al. (2006) result has a slight uptick in the 10–12 Gyr bin) and this may be a result of their samples being dominated by younger thin disk objects.

We also compare our SFH to that obtained by Rowell (2013) who implemented an algorithm to invert the WDLF in order to derive the local SFH. The resulting SFH, computed from the Harris et al. (2006) WDLF, is shown as the green-dotted–dashed line in the left-hand panel of Figure 10. Both SFHs are normalized, and given the varying units obtained by each study, only a comparison of the shape is appropriate. Both SFHs show two stages of star formation, however, we do not see a recent peak 2–3 Gyr ago, and the peak in our SFR attributed to the thick disk happens about 2 Gyr earlier. This is likely due to the larger volume sampled by our CFIS-PS1-*Gaia* sample, which in turn will allow us to better constrain this time frame because our sample will contain a larger fraction of thick disk objects. It is worth noting that similarly a burst in the SFR at 600 Myr is seen in the luminosity function of the local 40 pc white dwarf sample from Torres & García-Berro (2016), but not in the deeper luminosity function generated by Munn et al. (2017). Kilic et al. (2017) concluded that this discrepancy is a result of the local sample being confined to the midplane of the MW, where most recent star formation occurs.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Comparing three different realizations of the model with young ages (top row), intermediate ages (middle row), and old ages (bottom row) using color–color diagrams (left column) and RPMDs (right column). As the ages of each component increase, a larger number of white dwarfs at redder colors are visible due to the increased cooling ages, resulting in lower temperatures.</figcaption></figure>

In the right-hand panel of Figure 10, we compare our results to other SFHs found in the literature through GCE models. Toyouchi & Chiba (2018) used an open-box model, allowing for the inflow and outflow of gas, in order to reproduce the radial metallicity distribution obtained using APOGEE data in Hayden et al. (2015). This result was obtained using a sample of red-giant stars with $3 < R < 15$ kpc and $|z| < 2$ kpc, which is a much larger volume than our sample. Their result shows a rise in star formation beginning at 12 Gyr with a peak at 9 Gyr. This is followed by a roughly constant SFR to the present day. Our SFH matches the Toyouchi & Chiba (2018) results quite well at recent and ancient times but differs slightly at intermediate ages.

Also plotted are results by Snaith et al. (2015) who compared their GCE model to a sample of abundances acquired as part of the High Accuracy Radial Velocity Planet Searcher (HARPS) survey. They modeled the MW disk with two components: an inner ($R < 7$ kpc) closed-box and an outer ($R > 9$ kpc) component that can accrete gas. Their resulting inner disk SFH shows a two-phase formation scenario separated by a

period of inactivity lasting roughly 1 Gyr. A follow-up paper by Haywood et al. (2018) showed that the inner disk model is compatible with APOGEE data, and argued that the gap could be caused by the formation of a bar.

Their inner and outer disk results have been plotted in purple and orange, respectively, in Figure 10. Given the relative nature of their presented SFH, we have scaled their SFH so that both periods associated with the thin disk are equal. The inner disk shows an initial burst of star formation, which they attribute to the thick disk, and ends with a relatively constant SFR over the past 7 Gyr associated with the thin disk. Comparing the general shape of their inner disk SFH, an initial peak of SFH followed by a decline to a nearly constant SFR, shows that the two models are consistent, however, we do not *require* a gap in our SFH. This is likely because this is a feature of the inner disk, and because the gap is thought to be a result of the formation of a bar, it did not have a large effect on the local star formation we measure.

Our results differ in terms of the SFH of the thick disk, particularly at the earliest stages of the MW’s formation. This is
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Posterior distributions for the 13-dimensional parameter space sampled with *astroABC*. From left to right we show the mean functional age ($\xi$, Gyr), SFR ($M_{\odot}$ yr$^{-1}$), standard deviation ($\sigma_{\tau}$), and skew ($\alpha$) for the thin disk, thick disk, and halo. See Equation (3) for the definition of each parameter. The final histogram shows the fraction of white dwarfs with helium atmospheres, $f_{\text{He}}$.</figcaption></figure>

likely a consequence of the large uncertainties in both models at the earliest epochs of the MW and in particular the difficulty of detecting white dwarfs with cooling ages exceeding 11 Gyr. This could also be a result of our choice of scale height (See Section 5), as an increased scale height at these epochs would increase our SFR.

Haywood et al. (2019) argue that the solar neighborhood is better described by the outer disk evolution, given the enrichment history of the outer disk. In their outer disk scenario, their simulation includes a single accretion event at 10 Gyr (and thus the SFR between 10 and 13.8 Gyr is unconstrained) and results in a consistent SFR over the past

10 Gyr, similar to the inner disk. This is roughly consistent with our results, particularly within the past 8 Gyr.

### 4.3. Component Masses

The right-hand panel of Figure 9 shows the cumulative mass of the MW as a function of time, again with all 1000 particles shown in gray. This is done by integrating our SFH and density functions over the entire MW. We measure a total stellar mass of $(3.4 \pm 0.6) \times 10^{10} M_{\odot}$ for our thin disk, thick disk, and stellar halo. Also plotted is the result of Toyouchi & Chiba (2018), which is consistent within uncertainties with our result. Given
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Comparing our best-fit model (bottom panels) to the data (top panels) in both color–color (left column) and RPMD (right column). The best-fit ages are shown in the bottom-right panel.</figcaption></figure>

the relative nature of the other SFHs, we have chosen not to make a direct comparison with the results plotted in Figure 10.

The mass of the MW’s stellar disk is an important astronomical property, and as such, has been extensively studied. Our result is systematically lower than many previous results. Recent mass models of the MW have found a total disk mass between (3.6–5.5) × 10¹⁰ M☉ (see, e.g., Gerhard 2002; Flynn et al. 2006; McMillan 2011; Bovy & Rix 2013). It should be noted, however, that these estimates are highly dependent on the model assumptions, including the scale length and scale height, IMF, and whether or not a bulge/bar is included. Furthermore, it should be reiterated that our survey volume is rather small, and thus extrapolating our volume to the entire disk will inevitably contain a systematic uncertainty.

### 4.4. Local Stellar Density and White Dwarf Number Density

Given the systematic uncertainties associated with extrapolating our result to the entire disk, a more appropriate comparison is likely with the local stellar density. The stellar density as calculated by our model is (0.036 ± 0.004) M☉ pc⁻³

for the thin disk, (4.5 ± 0.5) × 10⁻³ M☉ pc⁻³ for the thick disk, and (3.2 ± 3.1) × 10⁻⁶ M☉ pc⁻³ for the halo. Our values agree with Bovy (2017), who used Gaia DR1 to calculate a main-sequence stellar density of 0.040 ± 0.002 M☉ pc⁻³.

The resulting white dwarf densities are (4.8 ± 0.4) × 10⁻³ pc⁻³ for the thin disk, (1.0 ± 0.2) × 10⁻³ pc⁻³ for the thick disk, and (6.3 ± 2.4) × 10⁻⁶ pc⁻³ for the halo. This corresponds to a total white dwarf number density of (5.8 ± 0.5) × 10⁻³ pc⁻³. Our values are consistent with the estimates of (5.5 ± 0.1) × 10⁻³ pc⁻³ from Munn et al. (2017), and marginally higher than results of 4.6 × 10⁻³ pc⁻³ from Harris et al. (2006), (4.8 ± 0.4) × 10⁻³ pc⁻³ from Torres et al. (2019), and (2.81 ± 0.52) × 10⁻³ pc⁻³ from Fantin et al. (2017).

Our resulting white dwarf contributions within the solar neighborhood break down to (83 ± 5)%, (17 ± 3)%, and (0.11 ± 0.05)% from the thin disk, thick disk, and stellar halo respectively. Rowell & Hambly (2011) constructed thin disk, thick disk, and halo luminosity functions using a statistical approach based on kinematics and photometry from the SuperCOSMOS survey in order to obtain local fractions of
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Left: positions of our model white dwarfs (red) within the MW (black). The concentric cyan rings represent heliocentric distances of 1, 2, and 3, kpc respectively. Right: heliocentric distance distributions of our model white dwarfs. Our sample consists mainly of local white dwarfs, with 96% having a distance less than 1 kpc.</figcaption></figure>

79%, 16%, and 5%, respectively. This result was consistent with Reid (2005) who found that thick disk white dwarfs should comprise roughly 20% of the local population. Our thick disk fraction is marginally lower than the recent result by Torres et al. (2019) who obtained fractional contributions of 74%, 25%, and 1% from the thin disk, thick disk, and halo, respectively.

### 4.5. Component Ages

The WDLF has been used as a tool to calculate the ages of various populations within the MW, including its field population. Our functional ages, $\xi$ (see Equation (3)), for the thin disk and thick disk are $8.41_{-0.35}^{+0.34}$ Gyr and $9.78_{-0.33}^{+0.36}$ Gyr, respectively. Because the age can be defined as the onset of star formation, we also present these values as $(8.5 \pm 0.3)$ Gyr for the thin disk and $(11.3 \pm 0.3)$ Gyr for the thick disk.

The first age measurement using white dwarfs was preformed by Winget et al. (1987) who obtained an age of $(9.3 \pm 2.0)$ Gyr for the Galactic disk. Following this, Leggett et al. (1998) used spectroscopic data from a sample of 43 white dwarfs to obtain an age of $(8 \pm 1.5)$ Gyr for the disk. Both of these samples were likely dominated by thin disk objects, and this is reflected in their age measurements.

Our SFH is consistent with those derived from the oldest open clusters, which are typically associated with the thin disk. Bedin et al. (2005) and García-Berro et al. (2010) used the WDLF to determine an age of ~8 Gyr for NGC 6791, an open cluster with solar-like metallicity. This cluster likely formed at the onset of thin disk star formation given its age and metallicity.

Our thick disk age is consistent with metal-rich globular clusters, like 47 Tucanae, which have metallicities comparable to thick disk field stars ($[Fe/H] = -0.75$). Hansen et al. (2013) used *Hubble Space Telescope* photometry to observe the WDLF and obtain an age of $(9.9 \pm 0.7)$ Gyr. In our SFH, this would place its formation near the peak of the thick disk star formation period.

Kilic et al. (2017) used the WDLF presented by Munn et al. (2017), who performed second epoch follow-up photometry within the SDSS footprint to calculate proper motions and select their white dwarf sample. The ages of our thin disk and

thick disk are consistent within the errors with their results of 7.4–8.2 Gyr and 9.5–9.9 Gyr.

Although less constrained than our disk ages, our resulting functional halo age is $10.92_{-1.09}^{+1.26}$ Gyr. This translates into an onset age of $(12.3 \pm 1.3)$ Gyr. This result is consistent with metal-poor globular clusters associated with the stellar halo, which typically have ages of 11–13 Gyr. For example, Hansen et al. (2013) obtained an age of $(11.7 \pm 0.3)$ Gyr for NGC 6397 and Bedin et al. (2009) found an age of $(11.6 \pm 0.6)$ Gyr for M4, both using the WDLF. Our result is also consistent with Kilic et al. (2017), who presented an age of $12.5_{-3.4}^{+1.4}$ Gyr.

### 4.6. He Fraction

The resulting He fraction in our sample is $(21 \pm 3)\%$. This is consistent with Bergeron et al. (2011), who fit model atmospheres to a spectroscopic white dwarf sample from the Palomar Green Survey. They showed that approximately 20% of white dwarfs with temperatures below 20,000 K appeared to contain a helium atmosphere. This result was followed up by Genest-Beaulieu & Bergeron (2019), who used SDSS spectra to show that the fraction of DB stars increases at lower temperatures (see their Figure 23). This fraction was found to be approximately 5% at temperatures greater than 20,000 K, with an increase to 20% –25% at 12,000 K. Using a sample of local white dwarfs Limoges et al. (2015) also showed that 25% of their sample was best described by models containing a helium atmosphere.

Our measurement is also consistent with recent results from Kilic et al. (2018) who used a CMD, generated using SDSS photometry and *Gaia* DR2 astrometry, to obtain a DB fraction of $36\% \pm 2\%$ for a sample of white dwarfs within 100 pc. Because He-atmosphere white dwarfs are most pronounced at higher temperatures, their result was quoted within a region of color–magnitude space encompassing the bifurcation of DA and DB white dwarfs due to the increased strength of the Balmer lines. In the temperature region, He-atmosphere white dwarfs are classified as type DB because He I lines are present in their spectra. Applying the same selection area, $M_{u_{SDSS}}$ between 10 and 14 and $u_{SDSS}-g_{SDSS}$ between –0.4 and +0.6, to our CFIS and PS1 model CMD produces a consistent DB fraction of $34\% \pm 3\%$.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Left: MW SFR as a function of lookback time. Right: the cumulative mass as a function of lookback time with the contribution from the thin disk (blue-dashed line), thick disk (green-dotted line) and halo (red-double-dotted–dashed line) highlighted.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

# 5. Discussion

As Figure 9 shows, our results suggest the star formation of the MW disk began (11.3 ± 0.3) Gyr ago. Following a brief period of halo formation, the thick disk rapidly begins to form stars until it reached a peak (9.8 ± 0.3) Gyr ago. This period is followed by a decline in thick disk star formation. The thin disk then began to form stars (8.5 ± 0.4) Gyr ago.

The transition between the thin and thick disk SFHs was a fairly active period in the MW's history. Recent results by Helmi et al. (2018) and Belokurov et al. (2018) provided evidence that this epoch was dominated by the accretion of at least one satellite galaxy, and may have contributed to the rapid transition between the thick and thin disk. The resulting merger would have heated up the existing stars to form a thick disk approximately 10 Gyr ago (Helmi et al. 2018). Given that a merger can enhance the SFR for a short period before quenching (Di Matteo et al. 2008), our results suggest that the merger occurred shortly before our peak SFR at (9.8 ± 0.3) Gyr, and before the onset of star formation in the thin disk. Thus, this scenario is consistent with a merger triggering the epoch of peak star formation in the MW and the start of star formation in what we now call the thin disk.

In this section, we discuss the impact of our assumptions on these results, including our choice of star formation prescription, and discuss the next steps for this study.

## 5.1. Effect of Star Formation Prescription

As shown in Equation (3), we assume skewed Gaussian functions for the three components, as allowing a skewness increases the degrees of freedom of the function. In this sense, we have assumed unimodal functions for the star formation histories, and we will inevitably fit over any potential multimodal peaks in the star formation. Specifically, if the SFH of a single component is composed of multiple bursts, then we will fit the average star formation throughout these periods of increased and decreased star formation. To examine this possibility, we have plotted the fractional residual between the data and the best-fit model in Figure 11.

Figure 11 shows a deficit of white dwarfs within the solid black box and an excess of white dwarfs within the dashed box.

These areas represent mean formation ages of (3.3 ± 1.8) Gyr and (5.8 ± 1.1) Gyr, respectively. This suggests that we are underestimating the SFR around 3 Gyr and overestimating at 6 Gyr by roughly 50% and 30%, respectively. This is consistent with the results of Rowell (2013) and Vergely et al. (2002), as shown in the left-hand panel of Figure 10, who saw an increase in star formation between 2 and 3 Gyr and a resulting decrease in the vicinity of 6 Gyr. This result is also in line with the recent work by Mor et al. (2019), who compared simulations of the Besançon Galactic model with the colors, magnitudes, and parallaxes from Gaia DR2, which shows a burst of star formation between 2 and 3 Gyr in the past. The slight bimodality suggested by Figure 11 is, however, not seen in the GCE models of Snaith et al. (2014) and Toyouchi & Chiba (2018), as seen in the right-hand panel of Figure 10.

## 5.2. Effect of the Scale Height

In our model, we have assumed a three-component Galaxy with constant scale heights. This assumption will inevitably impact our SFH, as a decrease in the scale height will result in a lower SFR. This is because more stars will reside closer to the plane, and will have a higher probability of passing our observational criteria.

Robin et al. (2014) found that the scale height for the thick disk (described using a sech² distribution) likely contracted from 800 to 350 pc between 12 and 10 Gyr. Given this scenario, our SFH would underpredict the SFR while the scale height was greater than 550 pc and overpredict while it was below this value. However, given the relatively small distances of the vast majority of our sample, the impact of deviating from our chosen mean of 550 pc is minimal. Specifically, at the maximum scale height of 800 pc we would need to produce 11% more stars than at a scale height of 550 pc. Similarly, at 350 pc we would produce 17% more stars than at 550 pc, reducing our SFR. These values are well within our uncertainties seen in Figure 9.

A thick disk composed of a range of scale heights has been suggested by recent results in order to explain the dip in [α/Fe] seen in the APOGEE data (Bovy et al. 2012; Haywood et al. 2013, 2016). Given that we modeled a single scale height, adding further components to our model would result in both scenarios
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Left: our SFH (black) compared to *Hipparcos* results by Vergely et al. (2002) (red-dashed line) and Cignoni et al. (2006) (cyan-dotted line). Also shown is the result of Rowell (2013), who inverted the WDLF of Harris et al. (2006) (green line). Right: comparing our SFH to results obtained via GCE models from Snaith et al. (2015) and Toyouchi & Chiba (2018). Due to the varying units presented by each study, the SFHs have been normalized, and therefore only the shape should be compared.</figcaption></figure>

<figure><img src="image_16.png" /></figure>

described above. The net result would be a flattening of the thick disk SFH with a minimal decrease to the peak because the total number of white dwarfs produced would need to remain the same.

### 5.3. Effect of Metallicity

In our model, we have chosen to set a constant metallicity for each Galactic component. The metallicity of a given star affects the pre-white dwarf lifetime through the analytical lifetimes from Hurley et al. (2000). Therefore, changing the metallicity of a given star will change its lifetime, and in turn its white dwarf cooling age. Specifically, for a fixed initial mass, the progenitor lifetime will increase with increasing metallicity, which in turn decreases the white dwarf cooling time. Tononi et al. (2019) found that including a dispersion in the metallicity provided a better fit to the 40 pc WDLF. Adding a dispersion to our metallicity would add a dispersion to the progenitor ages, which, through the cooling age, affects the photometry. If the dispersion is symmetric the effect would be equal at all ages and thus not impact our results. We have studied this effect in Figure 12, which shows the fractional change in progenitor lifetime for a higher and lower metallicity relative to our chosen values for the thin and thick disk. While not perfectly symmetric, the implementation of a constant dispersion would only systematically increase our ages by a few percent, which is minimal relative to our uncertainties.

### 5.4. Effect of the IMF and IFMR

Ultimately, as with any population synthesis study, our results are contingent on our assumptions. Assuming that the IMF is constant over time, changing our IMF to, for example, a Salpeter or Chabrier will affect the overall normalization, but not the shape of our SFH. This is because much of the discrepancy between IMFs is present at the low-mass end and the white dwarfs currently present in the MW were formed by stars with initial massed greater than 0.8–0.9 $M_\odot$, for which the slope of the IMF is well constrained.

Changing the IFMR will impact the ages of the white dwarfs, as the white dwarf cooling age is a strong function of the mass.

As the IFMR is generally estimated from empirical data, *Gaia* DR2 can be used to calculate the IFMR. Recent work by El-Badry et al. (2018), who used spectroscopically confirmed DAs in wide binaries, has shown promise, however, extending this work to all white dwarfs has been challenging given the impact of the atmospheric type on the resulting cooling age, which requires spectroscopic follow-up. The results by El-Badry et al. (2018) and Cummings et al. (2018) also showed that much of the discrepancy between IFMRs arises at high progenitor masses, for which there are very few such stars in our observed sample (see their Figure 3). While we note that there is no consensus on an IFMR, changing the IFMR is likely to have the highest impact on the SFH at old ages where the majority of old white dwarfs are formed from higher mass progenitors.

In order to test this, we ran a simulation with the semiempirical IFMR derived by Cummings et al. (2018) using the MESA Isochrones and Stellar Tracks (MIST) isochrones from Choi et al. (2016) to calculate the progenitor ages. The IFMR presented by Cummings et al. (2018) produces systematically larger white dwarf masses relative to the IFMR from Kalirai et al. (2008), with typical differences being ~0.05 $M_\odot$. At a fixed SFH, or equivalently, fixed white dwarf cooling time, increasing the white dwarf mass will increase the resulting temperature. The end result is that the ages must increase to match the observed colors as these white dwarfs will require more time to cool to equivalent temperatures, represented by the observed colors in our model. The resulting SFH with the Cummings et al. (2018) IFMR can be seen in Figure 13. The shape of the two SFHs are roughly consistent, although the resulting thick disk functional age is ~0.3 Gyr older than the result with the Kalirai et al. (2008) IFMR seen in Figure 9.

### 5.5. Effect of Atmospheric Composition

In our model, we assumed that white dwarfs have one of two atmospheric compositions: pure hydrogen or pure helium. In reality, white dwarfs with helium dominant atmospheres typically contain some hydrogen (Koester & Kepler 2015;
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Color–color diagram showing the data (left) and best-fit model (center) binned every 0.15 mags. The right-hand panel shows the difference between the data and the model, color-coded by either an excess (blue) or deficit (red) within the model. The black box shows the main location where the model contains a deficit of white dwarfs relative to the data, and these objects have a mean formation age of (3.3 ± 1.8) Gyr. The dashed box represents the location where our model overpredicts the number of white dwarfs, and these objects have a mean formation age of (5.8 ± 1.1) Gyr. This suggests a more bimodal formation history, with a 50% increase in SFR near 3 Gyr and a 30% deficit at 6 Gyr.</figcaption></figure>

Rolland et al. 2018). Genest-Beaulieu & Bergeron (2019) showed that this additional hydrogen impacts the photometric mass determinations of helium-atmosphere white dwarfs, particularly at low temperatures where the hydrogen abundance can increase due to convective mixing. Given that our model obtains masses from the IFMR (and hence both populations have the same mass distribution), this issue will not impact our results, however, a change in the atmospheric composition will affect the rate of cooling of the helium-atmosphere white dwarfs. As discussed in Kilic et al. (2017), pure-helium and pure-hydrogen white dwarfs take the same amount of time to cool to roughly 5000 K, however, hydrogen-atmosphere white dwarfs begin to cool slower than their helium counterparts below 5000 K. Specifically, for a $0.6 M_{\odot}$ white dwarf with a pure-hydrogen atmosphere, it will take ~1.5 Gyr longer to reach 4000 K than if it had a pure-helium atmosphere. At $0.8 M_{\odot}$, this discrepancy rises to ~2.5 Gyr. Within this temperature range, our resulting model produces very few white dwarfs, and those that pass our selection criteria have pure-hydrogen atmospheres. This is consistent with recent results suggesting that the majority of cool white dwarfs have nearly pure-hydrogen atmospheres (Kowalski & Saumon 2006; Giammichele et al. 2012; Limoges et al. 2015; Bergeron et al. 2019). The addition of mixed H/He atmospheres may allow a few He-dominant white dwarfs to pass the selection criteria by slowing down their cooling at this epoch, however, this effect would be minimal given the scarcity of objects in this region.

## 5.6. Model Improvements

A few straightforward improvements could be made to our white dwarf population synthesis model in the future. The first would be the inclusion of binaries within the model. Many main-sequence+white dwarf binaries would not pass our selection criteria given their colors, which would increase the resulting star formation because the model would need to produce more stars to match the data.

Binary evolution during previous stages of stellar evolution can also affect the evolution of white dwarfs. Recent evidence presented by Kilic et al. (2018) suggests that there exists a population of white dwarfs that formed through the merging of two main-sequence stars, and that this population may account

for as much as 15%–25% of the field white dwarf population. This result is consistent with the theoretical predictions made by Toonen et al. (2017) using the local 20 pc sample. Single white dwarfs formed through merging events have also been suggested by Liebert et al. (2005) to explain the high-mass bump in the DA white dwarf mass distribution, but at a rate of 12%–15%. This scenario will not only affect the progenitor age, but also the cooling age because higher mass white dwarfs cool more quickly than their lower-mass counterparts. This will inevitably affect the measured SFH, however, further work will be needed to determine the exact contribution from this scenario.

Extremely low-mass white dwarfs ($M < 0.3 M_{\odot}$) are also likely present within our data set. These white dwarfs cannot be formed through normal evolutionary channels and are thought to form as a result of extreme mass loss before the horizontal branch, leaving an exposed core composed primarily of helium (see, e.g., Kilic et al. 2011; Sun & Arras 2018). Kepler et al. (2007) found that this population could account for as much as 10% of the total white dwarf population, and they have been found in many white dwarf surveys (Kleinman et al. 2013; Kepler et al. 2016, 2019). Recent work by Pelisoli & Vos (2019) using Gaia DR2 suggests a space density of 275 kpc$^{-3}$, which would result in only a handful of objects within our sample.

## 5.7. Better Data: A Look Ahead

Finally, we reiterate the need for deeper proper motion surveys in order to constrain the properties of the Galactic halo population. Future surveys like the LSST, WFIRST, Euclid (Laureijs et al. 2011), and The Cosmological Advanced Survey Telescope for Optical and ultraviolet Research (CASTOR; Côte et al. 2012) will usher in a new era in deep wide-field surveys that will uncover thousands of halo white dwarf candidates.

In order to emphasize this point, we have run a simulation with the mean parameters acquired as part of our fit to the depth of CFIS *u* (*u* = 24.2). The resulting color–color diagram and RPMD can be seen in Figure 14, where the right-hand panel shows the simulation with the CFIS magnitude limit and the left-hand panel is the same as Figure 7. The resulting simulation with the CFIS magnitude limit returns half an order of magnitude more white dwarfs, and the disparity is particularly noticeable in the thick disk and halo populations. Specifically, we expect that our sample of
<figure><img src="image_18.png" /><figcaption>**Figure 12.** Fractional difference between the mean metallicity (black solid line) for the thin disk (left) and thick disk (right) for a star with higher metallicity (blue-dotted line) and lower metallicity (red-dashed line) than the mean.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 13.** A comparison between the resulting SFH using the IFMR from Kalirai et al. (2008) (dashed line) and Cummings et al. (2018) (solid line). The Cummings et al. (2018) IFMR returns systematically higher masses (~0.05 $M_\odot$), which results in a 0.3 Gyr increase in the age for the thick disk.</figcaption></figure>

halo objects would increase by a factor of 5–10, and our thick disk sample by a factor of 3 at the depth of the CFIS observations. Given their importance in understanding Galactic evolution, in particular the age of the halo via the WDLF, upcoming surveys such as LSST and WFIRST will provide an unprecedented look into the evolution of the Galactic inner halo via its white dwarf population.

## 6. Summary

In this paper, we have presented a newly developed white dwarf population synthesis code and used it to determine the

SFH of the Galactic thin disk, thick disk, and stellar halo using new exquisite u-band data from CFIS. Our population synthesis model takes MW geometry and extinction, as well as survey parameters such as completeness, into account in order to return a mock catalog of white dwarfs in the specified photometric bands. We use data from the CFIS, PS1, and Gaia DR2 to compare to our resulting simulations with a given SFH, which we have parameterized using skewed Gaussian functions. This sample was shown to be local to the solar neighborhood, with a median distance of 388 pc. Our main results are the following:

1. The SFH of the MW disk begins (11.3 ± 0.5) Gyr ago with the formation of the thick disk. The thick disk formed stars for 3 Gyr and reached a maximum of (8.8 ± 1.4) $M_\odot$ yr⁻¹ at (9.8 ± 0.3) Gyr.

2. The thick disk peak was followed by a decline in the star formation for 1 Gyr before the thin disk began forming stars at a roughly constant rate of (3.5 ± 0.3) $M_\odot$ yr⁻¹ for the past 8 Gyr. The maximum SFR of the thin disk occurred at (8.4 ± 0.3) Gyr.

3. Although the star formation shape parameters for the stellar halo are relatively poorly constrained, we find a mean age for the inner halo of 10.9±1.3 Gyr.

4. Studying the residuals reveals variations from a unimodal function, with an increase by as much as 50% at (3.3 ± 1.8) Gyr and a decrease of 30% at (5.8 ± 1.1) Gyr.

5. The resulting mass of the thin disk, thick disk, and stellar halo was found to be (3.4 ± 0.6) × 10¹⁰ $M_\odot$.

6. The white dwarf space densities were found to be (4.8 ± 0.4) × 10⁻³ pc⁻³ for the thin disk, (1.0 ± 0.2) × 10⁻³ pc⁻³ for the thick disk, and (6.3 ± 2.4) × 10⁻⁶ pc⁻³ for the halo. The fractional contributions of the thin disk, thick disk, and halo to the local white dwarf population were found to be (83 ± 5)%, (17 ± 3)%, and (0.11 ± 0.05)%, respectively.

7. We find a He fraction of (21 ± 3)%. This is consistent with previous studies of the solar neighborhood.
<figure><img src="image_20.png" /><figcaption>**Figure 14.** Comparing the resulting white dwarf populations based on the *Gaia* magnitude limit (left columns; *G* = 20.7) and the CFIS magnitude limit (right columns; *u* = 24.2). The fainter magnitude limit results in a (5–10)×increase in halo white dwarfs relative to the CFIS-PS1-*Gaia* sample.</figcaption></figure>

Given the *Gaia*-imposed magnitude limit of our survey, our data set is dominated by disk stars in the solar neighborhood. Future deep proper motion surveys will be needed to increase the sample of halo white dwarfs in order to properly constrain its formation parameters. Upcoming surveys, like Pan-STARRS1 DR2, the LSST, as well as *WFIRST*, *Euclid*, and *CASTOR* if they incorporate proper motions will open a new window into the formation of our MW at its earliest times.

We would like to thank the referee and all members of the CFIS for their insightful comments and contributions to this project.

This work was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC). ES gratefully acknowledges funding by the Emmy Noether program from the Deutsche Forschungsgemeinschaft (DFG). NFM gratefully acknowledges support from the French National Research Agency (ANR) funded project "Pristine" (ANR-18-CE31-0017) along with funding from CNRS/INSU through the Programme National Galaxies et Cosmologie.

This work is based on data obtained as part of the CFIS, a CFHT large program of the National Research Council of Canada and the French Centre National de la Recherche Scientifique. Based on observations obtained with MegaPrime/MegaCam, a joint project of the CFHT and CEA Saclay, at the CFHT that is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

The PS1 and the PS1 public science archives have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast,
the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

## ORCID iDs

Nicholas J. Fantin https://orcid.org/0000-0003-3816-3254

Patrick Côté https://orcid.org/0000-0003-1184-8114

Pierre Bergeron https://orcid.org/0000-0003-2368-345X

Rodrigo A. Ibata https://orcid.org/0000-0002-3292-9709

Raymond G. Carlberg https://orcid.org/0000-0002-7667-0081

Geraint F. Lewis https://orcid.org/0000-0003-3081-9319

Khyati Malhan https://orcid.org/0000-0002-8318-433X

Nicolas F. Martin https://orcid.org/0000-0002-1349-202X

Douglas Scott https://orcid.org/0000-0002-6878-9840

## References

Bedin, L. R., Salaris, M., Piotto, G., et al. 2005, ApJL, 624, L45

Bedin, L. R., Salaris, M., Piotto, G., et al. 2009, ApJ, 697, 965

Belokurov, V., Erkal, D., Evans, N. W., et al. 2018, MNRAS, 478, 611

Bergeron, P., Dufour, P., Fontaine, G., et al. 2019, ApJ, 876, 67

Bergeron, P., Wesemael, F., Dufour, P., et al. 2011, ApJ, 737, 28

Bovy, J. 2017, MNRAS, 470, 1360

Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115

Bovy, J., Rix, H.-W., & Hogg, D. W. 2012, ApJ, 751, 131

Brown, T. M., Smith, E., Ferguson, H. C., et al. 2006, ApJ, 652, 323

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Cignoni, M., Degl'Innocenti, S., Prada Moroni, P. G., et al. 2006, A&A, 459, 783

Côte, P., Scott, A., Balogh, M., et al. 2012, Proc. SPIE, 8442, 844215

Cummings, J. D., Kalirai, J. S., Tremblay, P.-E., et al. 2018, ApJ, 866, 21

Di Matteo, P., Bournaud, F., Martig, M., et al. 2008, A&A, 492, 31

Eggen, O. J. 1998, AJ, 115, 2397

El-Badry, K., Rix, H.-W., & Weisz, D. R. 2018, ApJL, 860, L17

ESA 1997, The Hipparcos and Tycho Catalogues, ESA SP-1200 (Noordwijk: ESA)

Fantin, N. J., Côté, P., Hanes, D. A., et al. 2017, ApJ, 843, 53

Farrow, D. J., Cole, S., Metcalfe, N., et al. 2014, MNRAS, 437, 748

Ferguson, A. M. N., Johnson, R. A., Faria, D. C., et al. 2005, ApJL, 622, L109

Flynn, C., Holmberg, J., Portinari, L., et al. 2006, MNRAS, 372, 1149

Fontaine, G., Brassard, P., & Bergeron, P. 2001, PASP, 113, 409

Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

García-Berro, E., Torres, S., Althaus, L. G., et al. 2010, Natur, 465, 194

Genest-Beaulieu, C., & Bergeron, P. 2019, ApJ, 882, 106

Gentile Fusillo, N. P., Tremblay, P.-E., Gänsicke, B. T., et al. 2019, MNRAS, 482, 4570

Gerhard, O. 2002, SSRv, 100, 129

Giammichele, N., Bergeron, P., & Dufour, P. 2012, ApJS, 199, 29

Gilmore, G., & Reid, N. 1983, MNRAS, 202, 1025

Girardi, L., Groenewegen, M. A. T., Hatziminaoglou, E., et al. 2005, A&A, 436, 895

Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2015, ApJ, 810, 25

Green, G. M., Schlafly, E. F., Finkbeiner, D., et al. 2018, MNRAS, 478, 651

Hansen, B. M. S., Anderson, J., Brewer, J., et al. 2007, ApJ, 671, 380

Hansen, B. M. S., Kalirai, J. S., Anderson, J., et al. 2013, Natur, 500, 51

Harris, H. C., Munn, J. A., Kilic, M., et al. 2006, AJ, 131, 571

Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132

Haywood, M., Lehnert, M. D., Di Matteo, P., et al. 2016, A&A, 589, A66

Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2013, A&A, 560, A109

Haywood, M., Di Matteo, P., Lehnert, M., et al. 2018, A&A, 618, A78

Haywood, M., Snaith, O. N., Lehnert, M. D., et al. 2019, A&A, 625, A105

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85

Holberg, J. B., & Bergeron, P. 2006, AJ, 132, 1221

Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543

Ibata, R. A., McConnachie, A., Cuillardre, J.-C., et al. 2017, ApJ, 848, 128

Kepler, S. O., Kleinman, S. J., Nitta, A., et al. 2007, MNRAS, 375, 1315

Kepler, S. O., Pelisoli, I., Koester, D., et al. 2016, MNRAS, 455, 3413

Kepler, S. O., Pelisoli, I., Koester, D., et al. 2019, MNRAS, 486, 2169

Jennings, E., & Madigan, M. 2017, A&C, 19, 16

Johnson, D. R. H., & Soderblom, D. R. 1987, AJ, 93, 864

Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864

Kalirai, J. S., Hansen, B. M. S., Kelson, D. D., et al. 2008, ApJ, 676, 594

Karim, M. T., & Mamajek, E. E. 2017, MNRAS, 465, 472

Kilic, M., Brown, W. R., Allende Prieto, C., et al. 2011, ApJ, 727, 3

Kilic, M., Hambly, N. C., Bergeron, P., et al. 2018, MNRAS, 479, L113

Kilic, M., Leggett, S. K., Tremblay, P.-E., et al. 2010, ApJS, 190, 77

Kilic, M., Munn, J. A., Harris, H. C., et al. 2017, ApJ, 837, 162

Kleinman, S. J., Kepler, S. O., Koester, D., et al. 2013, ApJS, 204, 5

Koester, D., & Kepler, S. O. 2015, A&A, 583, A86

Kowalski, P. M., & Saumon, D. 2006, ApJL, 651, L137

Kroupa, P. 2001, MNRAS, 322, 231

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193

Leggett, S. K., Ruiz, M. T., & Bergeron, P. 1998, ApJ, 497, 294

Liebert, J., Bergeron, P., & Holberg, J. B. 2005, ApJS, 156, 47

Limgos, M.-M., Bergeron, P., & Lépine, S. 2015, ApJS, 219, 19

LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv:0912.0201

Luyten, W. J. 1922, LicOB, 336, 135

Marjoram, P., Molitor, J., Plagnol, V., et al. 2003, PNAS, 100, 15324

McMillan, P. J. 2011, MNRAS, 414, 2446

Mor, R., Robin, A. C., Figueras, F., et al. 2018, A&A, 620, A79

Mor, R., Robin, A. C., Figueras, F., et al. 2019, A&A, 624, L1

Munn, J. A., Harris, H. C., von Hippel, T., et al. 2017, AJ, 153, 10

Noh, H.-R., & Scalo, J. 1990, ApJ, 352, 605

Norris, J., Bessell, M. S., & Pickles, A. J. 1985, ApJS, 58, 463

Pelisoli, I., & Vos, J. 2019, arXiv:1907.03766

Peng, X., Du, C., Wu, Z., et al. 2013, MNRAS, 434, 3165

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, arXiv:1807.06209

Reid, I. N. 2005, ARA&A, 43, 247

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130

Richer, H. B., Anderson, J., Brewer, J., et al. 2006, Sci, 313, 936

Robin, A. C., Bienaymé, O., Fernández-Trincado, J. G., et al. 2017, A&A, 605, A1

Robin, A. C., Reylé, C., Derrière, S., et al. 2003, A&A, 409, 523

Robin, A. C., Reylé, C., Fliri, J., et al. 2014, A&A, 569, A13

Rolland, B., Bergeron, P., & Fontaine, G. 2018, ApJ, 857, 56

Rowell, N. 2013, MNRAS, 434, 1549

Rowell, N., & Hambly, N. C. 2011, MNRAS, 417, 93

Snaith, O., Haywood, M., Di Matteo, P., et al. 2015, A&A, 578, A87

Snaith, O. N., Haywood, M., Di Matteo, P., et al. 2014, ApJL, 781, L31

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv:1503.03757

Sun, M., & Arras, P. 2018, ApJ, 858, 14

Thomas, G. F., McConnachie, A. W., Ibata, R. A., et al. 2018, MNRAS, 481, 5223

Tononi, J., Torres, S., García-Berro, E., et al. 2019, A&A, 628, A52

Toonen, S., Hollands, M., Gänsicke, B. T., et al. 2017, A&A, 602, A16

Torres, S., Cantero, C., Rebassa-Mansergas, A., et al. 2019, MNRAS, 485, 5573

Torres, S., & García-Berro, E. 2016, A&A, 588, A35

Toyouchi, D., & Chiba, M. 2018, ApJ, 855, 104

Tremblay, P.-E., Kalirai, J. S., Soderblom, D. R., et al. 2014, ApJ, 791, 92

van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJL, 771, L35

Vergely, J.-L., Köppen, J., Egret, D., et al. 2002, A&A, 390, 917

Winget, D. E., Hansen, C. J., Liebert, J., et al. 1987, ApJL, 315, L77