© ESO, 2014

1. Introduction

The Galactic halo is the oldest component of our Galaxy, containing metal-poor stars with high velocity dispersion. It contains a small percent of the total stellar mass of the Galaxy. Many questions about the formation of the halo and the Milky Way’s oldest stars are still to be answered, such as What is its star formation history (SFH)?, What is the initial mass function (IMF)?, and What fraction of halo stars resides in binaries? In this paper, we will investigate all three of these questions by studying the population of halo white dwarfs with a population synthesis approach.

White dwarfs (WDs) are an increasingly important tool used to study Galactic populations. Because they are the end product of low and intermediate mass stars, WDs are interesting objects of study for age determinations (eg. Hansen et al. 2007; Bedin et al. 2009). Since we have entered the era of large sky surveys, a huge amount of high-quality observational data of these stars is now or will soon become available. Physically, WDs are rather well understood, and they have been used as cosmic chronometers to study our Galaxy, as well as open and globular clusters, for more than two decades (Winget et al. 1987; see reviews by Fontaine et al. 2001; and Althaus et al. 2010, for example). Since in general, halo WDs are cool and faint, we confine ourselves to studying the ones in the solar neighbourhood.

The halo WD luminosity function (WDLF) was first derived by Liebert et al. (1989), based on six observed WDs with tangential velocities v_t exceeding 250 km s^-1. The most recent estimate is based on observations of 93 WDs with v_t> 200 km s^-1 in the SuperCOSMOS Sky Survey (Rowell & Hambly 2011, hereafter RH11). Theoretical halo WDLFs have been determined by, amongst others, Adams & Laughlin (1996); Isern et al. (1998); Camacho et al. (2007). Predictions for Gaia’s performances on WDs have been made by Torres et al. (2005). For a recent paper on this topic, see Carrasco et al. (2014). However, the effect of binary stars has never been studied in great detail. Furthermore, different initial parameters, stellar evolution codes, and WD cooling models were used in most of these papers.

For different assumptions about the IMF, SFH, and binary fraction, as well as for two different WD cooling models, we determine the WDLF and compare it with the observed halo WDLF in RH11. We derive both its shape and its normalization from an independent mass density of low-mass halo stars. We include not only single stars, but also focus on the contribution from WDs in binaries and WDs that are the result of a binary merger. Furthermore, we predict the properties of the population of binary WDs in the halo for a standard model and derive the number of halo WDs that can be detected by the Gaia satellite.

The setup of this paper is as follows: in Sect. 2 we explain our methods, in Sect. 3 we discuss our results, and our conclusions can be found in Sect. 4. In the concluding section we try to answer the question: What can Gaia observations of halo WDs teach us about the IMF, SFH, and binary fraction in the halo?

2. Model ingredients

We aim to derive the WDLF from first principles, i.e. not normalizing it to the observed WDLF, but using an independent estimate of the local stellar halo mass density to deduce a WDLF. A very important ingredient of our model is therefore the relation between this local density ρ₀, the stellar halo mass in the solar neighbourhood (the region that we simulate) and the IMF. In the next subsection, the expected number of halo stars is derived for three different IMFs. More details on this calculation can be found in the two appendices of this paper.

2.1. Initial Mass Functions

As a standard assumption, the IMF φ(m) can be written as a power law $\mathit{\phi }\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{\equiv }\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{m}}\mathrm{\propto }{\mathit{m}}^{\mathrm{-}\mathrm{(}\mathit{\gamma }\mathrm{+}\mathrm{1}\mathrm{)}}$(1)with N being the number of stars formed in the mass range m,m + dm and γ the IMF slope. We assume φ(m) to be independent of Galactic age or metallicity. Unless specified otherwise, N here represents the number of stars in the case that all stars are single (a binary fraction of 0). In Sect. 2.2 we explain how these numbers change with a nonzero binary fraction.

In a classical paper, Salpeter (1955) estimated γ = 1.35, and the corresponding IMF is nowadays referred to as a Salpeter IMF. Although not our standard model, one of the IMFs that we investigate in this paper is a Salpeter IMF for the whole mass range of stars (0.1 − 100 M_⊙). It is nowadays generally believed that the IMF flattens below 1.0  M_⊙, so this can be considered a bottom-heavy IMF.

In our standard model we split the IMF up into three power laws, following Kroupa et al. (1993): $\mathit{\phi }\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{\propto }\{\begin{array}{c}\\ \mathrm{35}\mathit{/}\mathrm{19}{\mathit{m}}^{-1.3}& \mathrm{if}\mathrm{0.1}\mathrm{\le }\mathit{m}\mathit{<}\mathrm{0.5}\mathit{,}\\ {\mathit{m}}^{-2.2}& \mathrm{if}\mathrm{0.5}\mathrm{\le }\mathit{m}\mathit{<}\mathrm{1.0}\mathit{,}\\ {\mathit{m}}^{-2.7}& \mathrm{if}\mathrm{1.0}\mathrm{\le }\mathit{m}\mathit{<}\mathrm{100.}\end{array}$(2)The thrid IMF that we investigate in this paper is the top-heavy IMF suggested by Suda et al. (2013). These authors argued that the IMF for stars with [Fe/H] < −2 is lognormal $\mathit{\phi }\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{\mathit{m}}\exp \left[\mathrm{-}\frac{{\log }_{\mathrm{10}}^{\mathrm{2}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]$(3)with median mass μ = 10 and dispersion σ = 0.4. Originally, this IMF was proposed by Komiya et al. (2007) for stars with [Fe/H] < −2.5 to explain the observed features of carbon enhanced metal poor stars, therefore we refer to it as the Komiya IMF. The higher metallicity stars would be formed according to a Salpeter IMF. Following the metallicity distribution function (MDF) of a two-component halo model (An et al. 2013), 24% of the zero-age main-sequence (ZAMS) stars with masses between 0.65 and 0.75  M_⊙ are formed according to a Komiya IMF. Therefore, when normalizing the WDLF properly, we expect more signatures from high-mass WDs (which cool fast and are thus faint) when choosing this IMF.

In order to determine the actual number of stars in a population N, one has to integrate φ(m), thereby setting the integration boundaries and the normalization constant. For example, integrating Eq. (1) with normalization constant A yields $\mathit{N}\mathrm{=}{\mathrm{\int }}_{{\mathit{m}}_{\mathrm{low}}}^{{\mathit{m}}_{\mathrm{high}}}\mathit{A}{\mathit{m}}^{\mathrm{-}\mathrm{(}\mathit{\gamma }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}\mathrm{)}}\mathrm{d}\mathit{m}\mathit{.}$(4)Hereafter, m_low = 0.1 and m_high = 100 for single stars and binary primaries. The value of A can be determined from an observed mass or number density of stars. We use the estimated local stellar halo mass density ρ₀ = 1.5 × 10^-4 M_⊙ pc^-3, based on observations of 16 halo stars in the mass range 0.1 ≤ m < 0.8 by Fuchs & Jahreiß (1998). These authors derived ρ₀ = 1 × 10^-4 M_⊙ pc^-3 as a firm lower limit. For a discussion on the correctness of this value compared to for example the lower estimate ρ₀ = 6.4 × 10^-5 M_⊙ pc^-3 (Gould et al. 1998), see Digby et al. (2003) and Helmi (2008). Since 0.8  M_⊙ is roughly the mass below which all stars can be considered unevolved, and 0.1  M_⊙ is our assumed lower mass boundary of all stars that are formed, this mass density is directly related to the total mass in unevolved stars M_unev in our simulation box. Our top-heavy IMF has two normalization constants, one for the very metal-poor stars and one for the higher metallicity stars. The normalization constants are derived in Appendix B.

Our simulation box represents the stellar halo in the solar neighbourhood, which we parameterize in a principal axis Cartesian coordinate system as (Helmi 2008) $\mathit{\rho }\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{0}}^{\mathit{n}}}{\left({\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{q}}^{\mathrm{2}}}\right)}^{\mathit{n}\mathit{/}\mathrm{2}}\mathit{,}$(5)with r₀ the distance from the Sun to the Galactic centre, q the minor-to-major axis ratio and n the power law exponent of the density profile. Throughout this paper, an oblate stellar halo (q < 1) is assumed. A sphere with radius ξ <r₀ around the Sun defines the minimum width and height of our simulation box. We show in Sect. 3.3 that ξ = 2.95 kpc is sufficient for our study of halo WDs. Furthermore, we choose r₀ = 8.0 kpc (Moni Bidin et al. 2012, an average of 16 literature measurements), n = −2.8 and q = 0.64 (Jurić et al. 2008). The simulated area with these values of r₀, ξ, n and q is shown in Fig. 1. We note that although n = −2.8 and q = 0.64 are the formal best-fit parameters of Jurić et al. (2008), one should keep in mind the ranges −3 ≤ n ≤ −2.5 and 0.5 ≤ q ≤ 0.8 as their fit results. Substituting the above-mentioned value of ρ₀ into Eq. (5) and integrating over the volume of our simulation box, we find M_unev = 3.6 × 10⁷ M_⊙ (see Appendix A).

Fig. 1 Simulation box containing a sphere with radius ξ = 2.95 kpc centred at the position of the Sun (x,y,z)⊙ = (8.0,0,0) kpc, n = − 2.8 and q = 0.64 (left panel) and a density map of the simulation box (right panel).

A crucial part of the normalization is the mass function of low-mass stars. The above-mentioned three IMFs predict drastically different numbers of stars in the range of masses 0.1 ≤ m < 0.8, which may or may not be in agreement with the observed sample from which the local halo mass density is determined (Fuchs & Jahreiß 1998). Since the mass function cannot be determined indisputably from this sample of 16 stars, we assume for each individual modelled IMF that it holds also in the low-mass regime. However, the resulting number of expected halo WDs depends strongly on the normalization of these low-mass stars (as we will see in Sect. 3.1), so we also investigate the effect of different slopes of the IMF at 0.1 ≤ m < 0.8.

Since the mass in low-mass stars is fixed by our normalization, the number of evolved stars and thus of WDs depends on the ratio of evolved stars to unevolved stars. The flatter the slope of the IMF for unevolved stars, the fewer unevolved stars there are, i.e. the higher this ratio¹. Most studies of low-mass stars suggest that the slope of this part of the IMF flattens (eg. Bonnell et al. 2007, γ_unev ≈ 0). This is why we take the Kroupa mass function as our standard. Furthermore, we calculate the number of evolved stars, which have initial masses 0.8 ≤ m < 100, for a flat (γ_unev = − 1) IMF, yielding a robust upper limit on the number of evolved stars, N_ev,upper (see Appendix B). In this way we derive a range of possible values for the number of evolved stars in our simulation box, between N_ev (where γ_unev is set consistently by the IMF, also in the low-mass regime) and N_ev,upper (derived by setting γ_unev = − 1). These numbers for the three different assumptions about the IMF are given in the first three columns of Table 1, assuming all stars are single. In Sect. 2.2 we give a discription of the last three columns of Table 1.

2.2. Population synthesis

The evolution of the halo stars is calculated with the binary population synthesis code SeBa (Portegies Zwart & Verbunt 1996; Toonen et al. 2012; Toonen & Nelemans 2013). In SeBa, stars are generated with a Monte-Carlo method, using the following distributions:

• Binary primaries are drawn from the same IMF as single stars (see Sect. 2.1).

• Flat mass ratio distribution between 0 and 1, thus for secondaries m_low = 0 and m_high = m_primary.

• Initial separation: flat in log a (Öpik’s law) between 1 R_⊙ and 10⁶R_⊙ (Abt 1983), provided that the stars do not fill their Roche lobe.

• Initial eccentricity: chosen from the thermal distribution Ξ(e) = 2e between 0 and 1 as proposed by Heggie (1975) and Duquennoy & Mayor (1991).

All simulated halo stars have metallicity Z = 10^-3 (0.05 Z_⊙), unless a top-heavy IMF is assumed. In that case, all stars that are born following a Komiya IMF are generated with metallicity Z = 2 × 10^-4 (0.01 Z_⊙). The common-envelope (CE) presciption of the standard model in SeBa (γα, Toonen et al. 2012) is used. With SeBa, we calculate the stellar evolution up to the point where the stars become WDs, neutron stars, or black holes, as well as the evolution of the binary systems until the end time of the simulation. For the WD cooling a separate code is used (see Sect. 2.4).

Having determined the total stellar mass in the simulated area, we still need to make an assumption on the binary fraction in order to arrive at the total number of stars in our simulation box. Because we assume a flat mass ratio distribution, the mass of the secondary is on average half the mass of the primary. The total number of binary systems in our simulation box if all stars are in binaries is therefore 1.5 times less than the total number of ZAMS stars if the binary fraction is zero. As a standard assumption we adopt a binary fraction of 0.5. This means that there are as many binary systems as there are single stars, thus that two out of every three stars are in a binary system. The total number of single stars (which is the same as the total number of binary systems) in this case can be found by dividing the numbers in the first three columns of Table 1 by 2.5.

The resulting number of WDs in our simulation box is listed in the last three columns of Table 1 for three different values for the binary fraction (0, 0.5, or 1) and standard assumptions about the SFH and WD cooling (see the next subsections). These numbers are quite sensitive to the assumed binary fraction, especially for a Kroupa or Salpeter IMF, because most of the binary primaries are unevolved stars in this case. This means that even more secondaries are unevolved stars, which do not become WDs within the age of the Galaxy, if the binary fraction is larger. Thus the total number of WDs is smaller if the binary fraction is larger. For a top-heavy IMF this effect is obviously less significant.

In our simulation we distinguish between carbon-oxygen core (CO) WDs, helium core (He) WDs and oxygen-neon core (ONe) WDs. He WDs must have undergone episodes of mass loss in close binary systems in order to be formed within a Hubble time. They are thus only found in binary systems, or as a result of two components of a binary system that merged. In Table 2, the mass ranges in which these three types of WDs occur in our simulation are listed. These mass ranges are partly overlapping, due to the effect of mass transfer in close binary systems. We note that they are dependent on the initial to final mass relation (IFMR), and therefore using a different population synthesis code may affect these results. See Toonen et al. (2014) for a comparison between four population synthesis codes.

Fig. 2 Comparison of cooling tracks of 0.2  M⊙ (dashed lines), 0.5  M⊙ (solid lines) and 1.2  M⊙ (dotted lines) WDs, between the models from the Althaus group (dark blue lines) and the Bergeron group (light blue lines). The 3 lines that represent cooling tracks from the Althaus models correspond to 3 WDs with different core compositions: He for the 0.2  M⊙ WD, CO for the 0.5  M⊙ WD and ONe for the 1.2  M⊙ WD, whereas the 3 lines from the Bergeron models all correspond to WDs with a CO core. The age of the universe is indicated with a vertical thin black line.

2.3. Star Formation Histories

We make three different assumptions about the SFH of the halo, based on observational indications that the vast majority of halo stars are old (Unavane et al. 1996; Kalirai 2012):

• (a)

  one burst of star formation 13 Gyr ago (ourstandard);

• (b)

  continous star formation from 13 until 10 Gyr ago, no star formation afterwards;

• (c)

  one burst of star formation 10 Gyr ago;

where a burst is assumed to last 250 Myr. After SeBa is run, all simulated halo stars have the same age, i.e. all stars are evolved for 10 Gyr or for 13 Gyr. To account for the SFH of the halo, we therefore shorten their lifetime by an amount of time randomly chosen between 0 and 250 Myr in case of assumptions (a) and (c) or with a number between 0 and 3 Gyr in case of assumption (b). If, by doing so, the lifetime of the star is reduced below the time it takes that star to become a WD, it is removed from our sample of WDs.

2.4. White dwarf cooling models

To determine the temperature, surface gravity and luminosity of a WD with a certain mass and cooling time, we use the cooling tracks published by Althaus et al. (2013) for He WDs and those from Renedo et al. (2010) for CO WDs. For ONe WDs, we use the cooling tracks from Althaus et al. (2007), both to determine their luminosities and temperatures, and their colours and magnitudes. We will refer to this set of cooling models as the Althaus models (our standard for WD cooling). The metallicities of the He and CO WDs in the Althaus models are assumed to be Z = 0.01, that of the ONe WDs Z = 0.02. Colours and magnitudes for CO WDs come from Althaus (priv. comm.), whereas colour tables of He WDs with high-metallicity progenitors (Z = 0.03) (Althaus et al. 2009) were used to determine the colours and magnitudes of He WDs. In all these cooling tracks and colour tables, the WDs have a higher metallicity than the ones in our simulation box (Z = 0.001). However, from the cooling tracks that were available for different metallicities (Althaus et al. 2009; Panei et al. 2007) we conclude that the effect of metallicity on WD cooling is smaller than other effects, such as the core composition (He or CO) of the WD, at least for large cooling times.

Alternatively, we also use the WD cooling tracks that are published online² (Bergeron et al. 2011; Holberg & Bergeron 2006; Kowalski & Saumon 2006; Tremblay et al. 2011), hereafter called the Bergeron models. The main difference between these two sets of cooling models is that in the Althaus models the evolution prior to WD formation is taken into account to arrive at a WD with a certain core composition, whereas in the Bergeron models the ad hoc assumption is made that all WDs have a CO core.

Fig. 3 From left to right: constructed grids of He, CO and ONe WD cooling in the Althaus models. The width of the panels corresponds to the maximum possible mass range per panel in our simulation, see Table 2. Vertical (horizontal) dashed lines indicate masses (cooling times) beyond which luminosities are obtained from extrapolation of the used cooling tracks. Dotted lines indicate the boundaries beyond which MV and MI magnitudes are obtained by extrapolation. For ONe WDs these regions completely overlap, for CO WDs the lower mass boundary partly overlaps, which is indicated by the light-grey-dark-grey dashed line. Scattered points indicate halo WDs that are observable with Gaia in our standard model. An explanation of the symbols is given Sect. 3.3.

In order to compare the two cooling models, we take the low-mass end of the CO WDs in the Bergeron models as the “He core” WDs in their models, and the high-mass end as their “ONe core” WDs, since non-CO core WD cooling models from Bergeron et al. were unfortunately not available in the literature. A comparison between these cooling models is given in Fig. 2. For details about the CNO flashes, which are very prominent on the cooling branch of the He WD in the left panel of Fig. 2, we refer to Althaus et al. (2013). Because the global specific heat of the He WDs is larger in the Althaus models, at a given cooling time the luminosity of such a low-mass WD is higher than that of a WD with the same mass in the Bergeron models (compare the two dashed lines in Fig. 2). Similarly, the global specific heat of ONe WDs is smaller in the Althaus models, resulting in lower luminosity WDs at a given cooling time compared to the high-mass WDs in the Bergeron models (the dotted lines in Fig. 2).

The explanation of the difference between the two solid lines in Fig. 2 is a bit more complicated, since in this case the WD core composition is the same in both models. Whether the prior evolution of the WD is or is not taken into account, will affect the onset of crystallization and the magnitude of the energy released by CO phase separation, a process that affects the cooling times below log (L/L_⊙) ≈ − 4. This could account for part of the difference in the cooling times between both models for the CO WDs. In addition, at low luminosities, the WD evolution is sensitive to the treatment of the outer boundary conditions and the equation of state at low densities (Althaus, priv. comm.).

Having chosen a set of cooling tracks, we want to determine the present-day luminosites for all the WDs in our model. Per WD type (He, CO, or ONe), typically only ten cooling tracks are available; these are interpolated and extrapolated over mass and cooling time to cover the whole parameter space that is sampled by our population synthesis code. The interpolation is linear, both in mass and in cooling time. For WDs that are more massive than the most massive WD of which a cooling track is available in the literature, we assume the cooling to be the same as for the most massive that is available. Similarly, the cooling track of the least massive WD is taken for WDs with a lower mass. In the Althaus models, this “extrapolation in mass” is done for He WDs with m < 0.155 or m> 0.435, for CO WDs with m < 0.505³ or m> 0.934, and for ONe WDs with m> 1.28. At the faint end of the cooling track, for CO WDs with m> 0.878 and t_cool ≳ 10¹⁰ years and ONe WDs with t_cool ≳ 6 × 10⁹ years (± 2 × 10⁹ years, depending on the WD mass), the luminosity is extrapolated using Mestel (1952) cooling. If the cooling tracks that we use in this study do not give a value for the luminosity of the WD at birth, we keep the luminosity constant at the value corresponding to the first given cooling time (for ONe WDs, this is ~10⁵ years). This yields lower limits to the luminosities of very young ONe WDs in the Althaus models and for all types of very young WDs in the Bergeron models.

In this way we construct three grids of WD cooling (one for each WD type; He, CO, or ONe), which are shown in the three different panels of Fig. 3. In this catalogue, the luminosity of every star in our simulation box can be found. The mass, type and cooling time of every WD in our simulation box was matched to the nearest catalogue grid point using the K3Match software (Schellart 2013). The dashed lines in Fig. 3 indicate the boundaries beyond which extrapolation was done as described above. Dotted lines similarly indicate the region beyond which the M_V and M_I magnitudes were determined using extrapolation. The luminosity and colour regions that are covered by the available cooling tracks in the literature overlap for ONe WDs, and so does the m = 0.505 boundary for CO WDs. This is indicated by the light grey dashes in between the dark grey vertical lines at the overlapping points. The scattered points in Fig. 3 visualize the position of the halo WDs that are observable with Gaia according to our standard stellar halo model. They will be discussed in Sect. 3.3.

2.5. Preparation of the WDLF

In this paper, we distinguish between spatially resolved and unresolved binaries. For each binary system, from the assigned distance and orbital separation a separation on the sky can be determined. Assuming thus two stars in a binary should be separated by at least 0.1−0.2 arcsec in order to be spatially resolved by Gaia (Arenou et al. 2005), we assign all binaries with a separation larger than or equal to 0.3 arcsec to the group of resolved binaries, those with a smaller separation to the group of unresolved binaries. Unresolved binaries are included as a single WD in the WDLF with a luminosity equal to the sum of the luminosites of the individual WDs in the binary.

To obtain a WDLF that can be compared with the observed one by RH11, we transform the luminosities of the WDs in our simulation box to bolometric magnitudes (using M_{bol, ⊙} = 4.75). We divide the total magnitude range into 2 bins per magnitude, ending up with bins such as M_{bol, ⊙} = [ 3.0,3.5 ] and [ 3.5,4.0 ], etc. The total number of stars per bin is then divided by the effective volume of our simulation box (V_eff = M_unev/ρ₀) to arrive at $\mathit{N}{\mathrm{pc}}^{-3}{\mathit{M}}_{\mathrm{bol}}^{-1}$.

There are also observational selection effects that need to be taken into account. Because RH11 only included halo WDs with tangential velocities v_t> 200 km s^-1, we reduce the number of WDs in each luminosity bin by a factor P(v_t> 200), which represents the probability that the tangential velocity of a halo star exceeds 200 km s^-1. The tangential velocities of halo stars (Chiba & Beers 2000) along the line of sight to one of the SuperCOSMOS survey fields is shown in Fig. 16 in RH11. From this figure, we estimate that P(v_t> 200) lies between 0.4 and 0.5, therefore we take 0.45. In our results section we will show the effect of choosing P(v_t> 200) = 0.4 or 0.5 instead.

2.6. Gaia magnitudes and extinction

The light from distant stars gets absorbed and reddened by interstellar dust. Following Toonen & Nelemans (2013), we assume that the dust follows the distribution $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{\propto }{\mathrm{sech}}^{\mathrm{2}}\mathrm{(}\mathit{z}\mathit{/}{\mathit{z}}_{\mathrm{h}}\mathrm{)}\mathit{,}& & \end{array}$(6)where z_h is the scale height of the Galactic dust (assumed to be 120 pc) and z the Cartesian coordinate in the z-direction. The interstellar extinction A_V between the Milky Way and a distant Galaxy in the V-band is assumed to be the extinction between us and a star at a distance d = ∞ (Sandage 1972), for Galactic latitude b = arcsin(z/d): $\begin{array}{ccc}{\mathit{A}}_{\mathit{V}}\mathrm{=}\left\{\begin{array}{c}\\ & \\ & \end{array}\right\}\mathrm{\equiv }{\mathit{A}}_{\mathit{V}}\mathrm{\left(}\mathrm{\infty }\mathrm{\right)}\mathit{.}& & \end{array}$The fraction of the extinction between us and a star at a distance d with Galactic latitude b ≠ 0 and this extinction is then $\begin{array}{ccc}\frac{{\mathit{A}}_{\mathit{V}}\mathrm{\left(}\mathit{d}\mathrm{\right)}}{{\mathit{A}}_{\mathit{V}}\mathrm{\left(}\mathrm{\infty }\mathrm{\right)}}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{d}\sin \mathit{b}}\mathit{P}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{d}\mathit{z}}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{P}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{d}\mathit{z}}\mathrm{=}\tanh \left(\frac{\mathit{d}\sin \mathit{b}}{{\mathit{z}}_{\mathrm{h}}}\right)\mathrm{·}& & \end{array}$(7)The stars in our simulation box are distributed according to the density profile given by Eq. (5), from which the distance to these stars is determined as $\begin{array}{ccc}\mathit{d}\mathrm{=}\sqrt{\mathrm{(}{\mathit{r}}_{\mathrm{0}}\mathrm{-}\mathit{x}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}\mathit{.}& & \end{array}$(8)Equation (7) is used to calculate the apparent magnitude V of a star at distance d and Galactic latitude b ≠ 0: $\begin{array}{ccc}\mathit{V}\mathrm{=}{\mathit{M}}_{\mathit{V}}\mathrm{+}\mathrm{5}\left({\log }_{\mathrm{10}}\mathrm{\left(}\mathit{d}\mathrm{\right)}\mathrm{-}\mathrm{1}\right)\mathrm{+}{\mathit{A}}_{\mathit{V}}\mathrm{\left(}\mathrm{\infty }\mathrm{\right)}\tanh \left(\frac{\mathit{d}\sin \mathit{b}}{{\mathit{z}}_{\mathrm{h}}}\right)\mathrm{·}& & \end{array}$(9)Because A_I(d) = 0.6 A_V(d) (Schlegel et al. 1998), we similarly calculate I from M_I and Eq. (7), after which Gaia magnitudes are calculated using (Jordi et al. 2010) $\begin{array}{ccc}\mathit{G}\mathrm{-}\mathit{V}\mathrm{=}{\mathit{a}}_{\mathrm{1}}\mathrm{+}{\mathit{a}}_{\mathrm{2}}\mathrm{(}\mathit{V}\mathrm{-}\mathit{I}\mathrm{)}\mathrm{+}{\mathit{a}}_{\mathrm{3}}\mathrm{(}\mathit{V}\mathrm{-}\mathit{I}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}_{\mathrm{4}}\mathrm{(}\mathit{V}\mathrm{-}\mathit{I}{\mathrm{)}}^{\mathrm{3}}& & \end{array}$(10)with a₁ = − 0.0257, a₂ = −0.0924, a₃ = −0.1623 and a₄ = 0.0090. We expect Gaia to detect all WDs with G < 20 (Brown 2013).

3. Results

We start our results section with an analysis of the theoretically determined WDLF in our standard halo model and compare it with the observationally determined one by RH11. In the second part of this section, we will compare the WDLFs predicted by models that were introduced in Sect. 2 and discuss our findings. In Sect. 3.3, we examine the halo WD population in more detail, again for our standard model. We derive the number of halo WDs that will be detectable by the Gaia satellite, and also predict properties of the whole population of (binary) WDs in the halo.

3.1. Standard model WDLF: theory vs. data

Fig. 4 Build-up of the WDLF in our standard halo model (50% binaries). All WDs are included in the black solid line (the total WDLF). The lower solid line shows the contribution from unresolved binary WDs. The dashed, dot-dashed and dotted lines show the contributions from CO, He and ONe WDs respectively. The light blue line with error bars is the halo WDLF determined from the SuperCOSMOS Sky Survey (Rowell & Hambly 2011), shown for comparison.

In Sect. 2.4 we have seen that the cooling tracks for He core WDs are quite different from those with CO cores, which in turn differ from ONe core WD cooling tracks. The effect of this can be seen partly in Fig. 4, where we present how the WDLF for our standard halo model is built up from contributions of the various WD types. We first note the monotonic increase in the WDLF, which occurs because the cooling of WDs is a simple gravothermal process (Isern et al. 2013). The drop in the number of stars at M_bol ≈ 16 is a consequence of the finite age of the universe. As was suggested by e.g. Winget et al. (1987), the observation of this drop can be used to constrain the age of our Galaxy. A second peak in the WDLF (the dotted curve around M_bol ≈ 19 in Fig. 4) is expected to consist of ONe WDs, due to their fast cooling times, as was first pointed out by Isern et al. (1998).

The contribution from the He WDs to the WDLF is shown with a dot-dashed line in Fig. 4. These He WDs have an unseen neutron star (NS) or black hole (BH) companion or they are the resulting merger product of two stars in a binary. However, there are many more He WDs that contribute to the WDLF: those in unresolved binary WD pairs (the lower solid line in Fig. 4). Since they have slow cooling times, the contribution of He WDs to the WDLF is largest at the bright end. The unresolved binaries that end up in the second peak of the WDLF (around M_bol ≈ 18.5) are systems in which at least one of the two WDs has an ONe core. The main contributors to the WDLF are CO WDs, visualized with a dashed line in Fig. 4, which is just below the black solid line. These can be single CO WDs, CO WDs in wide binary WD pairs, but also CO WDs with a NS or BH companion or merger products. WDs with a main-sequence star as companion are not included in the WDLF, because the light from the main-sequence star will dominate the spectrum in that case.

Fig. 5 Comparison of halo WD luminosity functions corresponding to different assumptions about the IMF (top left), WD cooling (top right), binary fraction (bottom left) and SFH of the halo (bottom right). The WDLF corresponding to our standard model, indicated by the blue solid line in all panels, is constructed using the 1.2 × 107 WDs in our simulation box (see Table 1).

Figure 4 shows that our standard model WDLF lies below the observed WDLF (RH11; the light blue line with error bars in Fig. 4), however we shall see in the next subsection that this discrepancy disappears when we vary the normalization. Our integrated standard model luminosity function (the black solid line in Fig. 4) yields n_{Halo WDs} = 2.08 × 10^-5 pc^-3. This value is lower than the integrated value of the RH11 WDLF, (1.4 ± 5.6) × 10^-4 pc^-3, mainly because of their higher estimate of the number of WDs in the luminosity bins around M_bol ≈ 17. Our models predict that there are no WDs in these bins. Although the present-day estimate of the number of WDs with M_bol ≈ 17 should be regarded as an upper limit because of the large error bars, future observations on the shape of this faint end of the WDLF will help to constrain WD cooling models and the SFH of the halo, whereas the normalization of the WDLF, especially at the faint end, will help to constrain the IMF and binary fraction (see Sect. 3.2).

When comparing our theoretically determined WDLF with the observed one by RH11 apart from the missing faint end (which is not reached by SuperCOSMOS because it is a magnitude-limited survey), also the missing bright end catches the eye. Here, another selection effect plays an important role:

RH11 only included WDs with a tangential velocity larger than v_t,min = 200km s^-1, to filter out thin and thick disk WDs. Due to the mean lower proper motion completeness limit μ_min = 40mas yr^-1 across the SuperCOSMOS Sky Survey, the sample of RH11 is becoming less complete at a distance of approximately $\frac{{\mathit{d}}_{\max }}{\mathrm{pc}}\mathrm{=}{\left(\frac{{\mathit{p}}_{\min }}{\mathrm{arcsec}}\right)}^{-1}\mathrm{=}\frac{{\mathit{v}}_{\mathit{t,}\min }}{\mathrm{4.74}\mathrm{km}{\mathrm{s}}^{-1}}{\left(\frac{{\mathit{\mu }}_{\min }}{\mathrm{arcsec}{\mathrm{yr}}^{-1}}\right)}^{-1}\mathrm{\approx }\mathrm{1}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}}\mathrm{·}$(11)Here, p_min is the minimum parallax that is reached and we have used that a proper motion of 1 arcsec yr^-1 corresponds to a tangential velocity of 1 AU yr^-1 = 4.74 km s^-1 at 1 pc. Because at distances larger than ~1 kpc, young and bright halo WDs contribute more to the WDLF than fainter WDs, the bright end of the WDLF is not reached by SuperCOSMOS.

We expect this latter bias to be resolved by Gaia, which can do microarcsecond astrometry and, as we will show in Sect. 3.3, is expected to detect intrinsically bright WDs to distances of ~2.5 kpc. Although the bright end of the WDLF has already be determined from an empirical measure of the WD cooling rate in a globular cluster (Goldsbury et al. 2012), we will soon have a new window on the Galactic halo, when we start to explore bright field halo WDs with Gaia.

3.2. Comparing the WDLFs of different halo models

Eight different model WDLFs are visualized in Fig. 5, comparing the effect of a different IMF, cooling track, binary fraction and SFH. As explained at the beginning of Sect. 2, we calculate not only the shape of the WDLF, but also derive its normalization from an independent mass density estimate of halo stars. For each model in Fig. 5 a band is given rather than a single line, which comes from the two different normalizations explained in Sect. 2.1. To arrive at the upper lines, the lower line is simply multiplied with a normalization factor corresponding to the used IMF (5.9 / 1.9 for the Kroupa IMF, 6.7 / 1.1 for the Salpeter IMF and 330 / 80 for the top-heavy IMF, see Table 1). The blue band in each panel represents our standard model, labelled “Kroupa”, “Althaus”, “50% binaries, 50% singles” and “13 Gyr burst” respectively.

The left panels of Fig. 5 show that a different IMF or binary fraction affects the normalization of the WDLF, as expected from Sects. 2.1 and 2.2. Regarding the shape of these five different WDLFs in the left panels of Fig. 5, the differences are the largest at the extremely faint and the extremely bright end. We see that the WDLF corresponding to a model with a larger binary fraction resembles more closely the shape of the lower solid line in Fig. 4, where the contribution to the WDLF from unresolved binary WDs is shown.

From the top right panel in Fig. 5 it is clear that there is a significant difference in the shape of the faint end of the WDLF when a different assumption is made about WD cooling. The drop between the two peaks of the WDLF is less prominent when the Bergeron models are used compared to the Althaus models, because CO WDs with a luminosity log (L/L_⊙) < −4 cool faster in the Bergeron models than in the Althaus models (see the right panel of Fig. 2).

The logarithmic scale on the vertical axis implies that it will be observationally challenging to distinguish the three different models of SFH (shown in the bottom right panel). These differ slightly from each other at the faint end of the WDLF. As expected, the WDLF of a 10 Gyr halo drops off at lower magnitudes than a 13 Gyr old halo. Furthermore, the gap between the two peaks of the WDLF is more prominent in the models with a SF burst compared to models with a continuous SFH.

There is a quite good overall agreement between our theoretically predicted WDLFs and the observed one by RH11, except for the model with a top-heavy IMF, which overpredicts the number of WDs per luminosity bin at the faint end of the WDLF. This also follows from the reduced χ²-test that we conducted to compare the agreement between the different model WDLFs in Fig. 5 with the observationally determined WDLF quantitatively (see Table 3).

From the first column of Table 3 we see that our standard model and the model with 100% singles have the lowest reduced χ² values (χ² = 2.29 and 2.28 respectively), closely followed by the models with alternative SFHs (χ² = 2.31 for the model with uniform SF between 10 and 13 Gyr, and χ² = 2.35 for the model with a SF burst 10 Gyr ago). The fact that all these values are so close together can also be determined from Fig. 5, where these four curves almost completely overlap. The models with 100% binaries or a Salpeter IMF do slightly worse, due to their lower normalization. In all cases the line corresponding to the upper limit of the number of stars has worse agreement with the observed WDLF (${\mathit{\chi }}_{\mathrm{upper}}^{\mathrm{2}}$; second column) than the lower line corresponding to that same model. This seems to indicate that the low-mass part of the IMF does not turn over at ~1.0 M_⊙ to become completely flat, but rather has a negative slope.

We varied the normalization of the WDLFs by multiplying them with a free parameter f, to see how well we can fit the shape of the WDLF. We kept f as a free parameter, because there are many ways in which we could adapt the normalization, for example choosing a different γ_unev (as we did for calculating the upper limits), a different binary fraction or a different mass density in unevolved stars ρ₀. The results of this analysis (summarized in the parameter ${\mathit{\chi }}_{\min }^{\mathrm{2}}$) are given in the third column of Table 3 for each model with the corresponding f value in the fourth column. Without normalizing the model WDLFs, the model with 100% binaries comes out best, with a reduced ${\mathit{\chi }}_{\min }^{\mathrm{2}}$ value of 2.25. Although these minimum ${\mathit{\chi }}_{\min }^{\mathrm{2}}$ values lie close together for most of the models, the Top-heavy and Bergeron models still have the worst agreement with the WDLF observed by RH11. For the two models with alternative SFHs and the model with 100% singles the ${\mathit{\chi }}_{\min }^{\mathrm{2}}$ values are larger than the χ² value corresponding to our preferred normalization, because there is one degree of freedom less if we fix the normalization of the WDLF.

The χ² values are also affected by our assumption of P(v_t> 200). If we had chosen the value P(v_t> 200) = 0.4 or 0.5, our standard model χ² value would change to 2.39 or 2.22 respectively, and how this other choice affects the other curves can be determined from the parameter f. If f is larger than 1, the larger value P(v_t> 200) = 0.5 would reduce the χ² value, if f is smaller than one P(v_t> 200) = 0.4 would yield a better match.

Fig. 6 Properties of unresolved binary WD pairs with G < 20 for our standard model. Shown are double He WDs (circles), double CO WDs (squares), CO+He systems in which the He WD is the brightest (upward pointing triangles), He+CO systems in which the CO WD is the brightest (downward pointing traingles), pre-WDs (filled stars), and “other” stars, which are pre-WDs that are indistinguishable from MS stars or giants. After the label discriptions in the legend, the total number of WD binaries of that particular kind is given. Due to the low number of halo WDs with G < 20 we expect to find in the Milky Way, there is some statistical noise in this figure.

Fig. 7 Properties of all 1.3 × 106 unresolved binary WD pairs in our simulation box (maximum magnitude G = 35), for our standard model. Contour lines mark the regions in which 33%, 67% and 95% of the binaries are located. Higher density regions are given darker colours. The 1% binaries in the tail of the distribution are indicated by the scattered points. Compared to Fig. 6, this figure hardly has statistical noise because of the large number of halo WDs we expect to find in the Milky Way halo if we would go up to magnitude G = 35.

3.3. Halo white dwarfs detectable by Gaia

In this subsection of our results, we take a closer look at the population of halo WDs in our standard model and what fraction of this population can be seen by the Gaia satellite.

An important point to keep in mind when studying halo WDs, is that one is biased towards young and bright WDs in a magnitude-limited survey. Since the bright part of the WDLF is to a large extent built up by unresolved binary WD pairs (see Fig. 4), we first look at their properties. Figure 7 shows the properties of all unresolved binary WD pairs in our simulation box, whereas Fig. 6 focusses on the ~300 unresolved binary WD pairs with G < 20. We note that this also includes binaries with large orbital separations which have never undergone interaction, because at large distances these can still be unresolved. Of the two WDs in each binary, the properties of the brightest are plotted. A distinction is made between CO+He WDs and He+CO WDs, the second of the two WD types in each group is the brightest WD in the system. For most of the systems, this is also the youngest WD. However, in some CO+He systems the He WD was formed first and is still brighter than the later formed CO WD, which is possible since He WDs can be intrinsically brighter than CO WDs at birth and they in general have longer cooling times (see Sect. 2.4). In the legend of each panel in Fig. 6, the number of systems of that particlar kind is given in brackets. Due to the low number of halo WDs we expect to find in the Milky Way, there is some statistical noise in this figure. What we want to show here, are the global positions of the WDs in this diagram, without focussing on their individual positions.

A particular aspect of the cooling models that we use as our standard, is that the luminosity of He WDs stays constant for a long time (10⁵ −10⁹ yr, depending on the mass), before cooling starts. This can be seen from the dark dashed line in the left panel of Fig. 2. As a consequence of this feature, He WDs that are on this part the cooling track will be seen more often than CO WDs with the same cooling time. We will refer to them as pre-WDs, to indicate that these objects do not look like standard WDs, because they are brighter and have smaller surface gravities.

In the top panels of Figs. 6 and 7, which are log g − log T_eff diagrams, the pre-WDs (plotted with star symbols) are clearly distinguishable from the other WDs because of their low surface gravities. We define pre-WDs as those double WDs in which the brightest of the two has 4.3 < log g < 6. There are more He WDs with even lower surface gravities, indicated by the label “other” in Fig. 6. However, these will be hard to distinguish from main-sequence stars or giants, which lie on or above the dashed line at log g ≈ 4 in Fig. 6 (Allen 1973). Because pre-WDs are only apparent when we use the Althaus cooling tracks for He WDs (see Fig. 2), Figs. 6 and 7 would not have any points with log g < 6 if we use the models with Bergeron cooling instead of our standard model. In the top panel of Fig. 7 we see that more than 95% of the halo WDs have log g> 7.0 and log T_eff < 4.2, which is also the part of the diagram where most of the single WDs and resolved WD binaries are expected to be situated. Furthermore, we see from the top panel of Fig. 6 that there is a narrow gap between the log g and log T_eff values of systems in which the brightest WD has a CO core and those in which the brightest of the two has a He core. In this way, these systems can in principle be distinguished from their positions in the log g − log T_eff diagram.

The middle panels of Figs. 6 and 7 show the IFMR for halo WDs, i.e. the mass of the brightest WD in the binary system (M_bright) is plotted as a function of its corresponding initial mass (M_ZAMS). In Fig. 6 we see that in most of these unresolved binary WDs the brightest of the two stars has a main-sequence progenitor star with a mass M_ZAMS ≈ 0.84 M_⊙. Because in our standard model halo stars are born about 13 Gyr ago, these stars have just become WDs, thus will be very abundant in the Gaia catalogue of halo WDs. There are very few high-mass WDs in our sample, mainly because their progenitor stars have shorter main-sequence lifetimes and they have thus cooled down much more. In the middle panel of Fig. 7, the global IFMR for halo WDs is shown. There is no focus on only the brightest WDs, with the result that the highest density region is shifted towards a line that resembles the IFMR of single stars, populated by the unresolved binary WDs that have not undergone interaction.

In the bottom panels of Figs. 6 and 7 the M_bright − orbital period relation for halo WDs is shown. We see that the unresolved binary WD pairs with G < 20 lie on three distinct lines, where each line is mostly populated by one of the different binary types. The majority of double CO systems have not interacted and thus evolve to systems with wide periods. Because most stars are low-mass they form WDs with similar masses, while the periods are determined by the initial period distribution. The short-period branch shows systems that are formed via CE evolution. In our model, this CE between a giant and a WD is always described by a the energy balance (α, see Toonen et al. 2012). The correlation they show between WD mass and orbital period can then be understood from the relation between the core mass and the radius of giants. Systems that start the CE phase in a more compact orbit will have giants with smaller radii and thus lower-mass cores. This means both a spiral in to shorter final periods and a final WD mass that is lower. The branch with longer periods shows systems that are formed via a second phase of mass transfer that was stable. During the mass transfer the orbit widens, which stops when the whole envelope of the giant has been transferred to the first formed WD. Due to the same relation mentioned above, giants with larger core masses (that form more massive WDs) are bigger and thus end their evolution in binaries with longer orbital periods. The same relation is seen in the WD companions to millisecond radio pulsars (Savonije 1987). From Fig. 7 it is clear that Fig. 6 only resembles a small part of the complete parameter space, but it constitutes a representative selection of the low-mass part of this diagram.

White dwarfs in unresolved binaries are of course not the only halo WDs we expect Gaia to observe. As we already mentioned in Sect. 3.1, single WDs, resolved double WDs, WDs with a NS or BH companion, and WDs that are the result of a merger also contribute to the WDLF. The number of WDs in each of these five groups is specified per WD type (pre-WD, He, CO, or ONe core) in Table 4. We see that all single WDs and all brightest WDs in a resolved binary system have a CO core with the limiting magnitude G < 20. If we look at fainter magnitudes, e.g. G < 23 or G < 26, Table 4 shows that ONe WDs will be detected, although there are still very few of them compared to CO WDs. The same is true for WDs with a NS or BH companion.

When selecting halo WDs from the Gaia catalogue, selection effects are expected, like the factor P(v_t> 200) that we multiplied our theoretically determined WDLFs with in the previous subsections to compare our results with that of RH11. Here, we do not include this factor, since it is not yet clear how large it will be. For example, for some fraction of the stars (V < 17), radial velocities will also be available. Therefore it should be possible to obtain a larger number of halo WD stars than just with a cut in v_t. Furthermore, the determination of the initial number of stars in our simulation box has a greater effect on the number of halo WDs than these selection effects have.

In the top and bottom rows of Table 4, the total number of halo WDs in two spheres around the Sun with respective radii 400 pc and 2.95 kpc is given, as well as the total mass these halo WDs constitute. From these we calculate the number densities of halo WDs n_{Halo WDs} = 5.89 × 10^-5 pc^-3 (within 400 pc) and the slightly higher value n_{Halo WDs} = 6.00 × 10^-5 pc^-3 (within 2.95 kpc). These values are more than a factor of two lager than the number density we derived by integrating the WDLF in Sect. 3.1. This difference is due to the factor P(v_t> 200) which is not taken into account here. Furthermore, here all halo WDs are counted within spheres of a certain radius around the Sun, whereas in Sect. 3.1 we estimated the number density including the edges of our simulation box (which is not a sphere, see Fig. 1).

Fig. 8 Distances and bolometric magnitudes of the 1.5 × 103 halo WDs that can be observed with Gaia, for our standard model. Top panel: distribution of their distances. Right panel: distribution of their bolometric magnitudes, which gives an idea of the statistical errors that are to be expected per luminosity bin of the WDLF for Gaia. The yellow stars, labelled “other”, are pre-WDs that are indistinguishable from MS stars or giants. These are not included in the projected distribution histograms.

Torres et al. (2005) also estimated the number of (single) halo WDs with G < 20 within 400 pc, and found 542 (their Table 3). We found slightly more halo WDs within 400 pc: 621, including both singles and binaries, see the top row of Table 4. However, it is not strange that these numbers differ from each other, given the large number of uncertainties in our estimate of the number density of halo WDs from the observed mass density in unevolved stars (Sect. 2.1) and the selection effects. The latter are implicitly taken into account by Torres et al. (1998), because they normalized the number density of halo WDs within 400 pc to the local observed value (Torres et al. 1998).

To check that our simulation box is large enough and make the claim that Gaia can detect approximately 1.5 × 10³ halo WDs with G < 20 (Table 4), we plotted the distances to all these WDs as a function of their bolometric magnitude in Fig. 8. With the same markers as in Fig. 6 the unresolved binaries are visualized, additional markers are used for single WDs, resolved binary WDs and WDs that are merger products. We see that apart from a few outlying “other” WDs, all WDs fit well within the sphere with radius ξ = 2.95 kpc around the Sun, which validates the size of our simulation box. As explained above, the “other” WDs were excluded from our luminosity function anyway because of their low surface gravities.

Projected onto the vertical axis of Fig. 8 is the number of halo WDs that can be found in every luminosity bin of the WDLF from the Gaia catalogue. From this right panel of Fig. 8 we get an idea of the statistical errors that are to be expected per luminosity bin of the WDLF. We see that the faint end of the WDLF will stay underdetermined since we expect to detect only a handful of WDs with M_bol> 16 with Gaia in our standard halo model. However, already with the few WDs in the lowest luminosity bins that can be reached, we can start comparing our halo models. It is not clear whether the drop at M_bol ≈ 16 is a detection limit, since a cut-off of the luminosity function due to the age of the Galaxy is expected at approximately the same bolometric magnitude (see Fig. 4).

As we explained at the end of Sect. 3.1, most of the WDs at the bright end of the WDLF can be included with Gaia whereas they could not before, since Gaia has a lower mean lower proper motion completeness limit than previous surveys. From the long tail of the distance distribution (top panel of Fig. 8), we see that there are many halo WDs with G < 20 beyond ~1 kpc, which all have absolute bolometric magnitudes M_bol < 8. It is because of the inclusion of these WDs that the bright end of the WDLF will probably be better constrained with Gaia than ever before.

The masses and cooling times of the 1.5 × 10³ halo WDs with G < 20 are plotted on top of the interpolated cooling track panels of Fig. 3. Again, the same markers are used as in Figs. 6 and 8. For the halo WDs that are merger products, it is indicated what type of WD the merger product is (He WD or CO WD) by the particular panel the marker is drawn in. Plusses and diamonds represent single WDs and resolved binary WDs respectively, which all have a CO core, as can also be seen from Table 4. This is not surprizing since He WDs can only be formed through binary interaction within the age of the universe and ONe WDs cool so fast that they pile up at the faint end of the WDLF, which will not be covered by Gaia.

It is interesting to see the narrow line at m = 0.54 M_⊙ in Fig. 3, where the single and non-interacting binary WDs pile up. This can be explained by the evolution lifetime of single and non-interacting ZAMS stars with an initial mass of 0.84 M_⊙, which is equal to the age of the halo in our standard model and we see in the middle panel of Fig. 7 that they become 0.54 M_⊙ WDs. Kalirai (2012, hereafter K12) first pointed out that the mass determination of these bright single halo WDs can be used to determine the age of the inner halo. The determination of the masses of the brightest WDs in a globular cluster provides an anchor point on the IFMR for low metallicity stars, since their ages can be deduced independently from the cluster age. K12 then drew a straight line through this anchor point and the mass of the brightest halo WDs, yielding a linear IFMR for WD masses between 0.50 and 0.58 M_⊙, see Fig. 9. The age of the field halo stars was subsequently determined using the Dartmouth Stellar Evolution Database (Dotter et al. 2008), in which ~0.55 M_⊙ WDs (with 0.83 M_⊙ progenitors) are ~11.4 Gyr old.

Fig. 9 Age of the brightest single halo WDs (and thus of the Galactic halo) as a function of WD mass. Since for single stars every age corresponds one-to-one to an initial stellar mass and WD mass, this is an alternative representation of the IFMR. The IFMRs for single stars predicted by SeBa (light blue curves) and the those predicted by MESA (dark blue curves) both differ from the black solid straight line through the data points with error bars of Kalirai (2012). The dashed lines correspond to the metallicity value that we used as our standard for the Milky Way halo in this study (Z = 0.001), while the dotted lines correspond to half that metallicity value (Z = 0.0005).

Although K12 did a carful analysis on the binarity of field halo WDs and in the globular cluster, the binary fraction in globular clusters is generally lower than in the field (eg. Ji & Bregman 2013, and references therein). Especially for the latter, the observed single WDs could thus still have an unseen companion or be the result of a binary merger. Furthermore, the behaviour of the (metallicity-dependent) IFMR in this mass regime is not well determined yet theoretically. This can be seen from Fig. 9, where we compare the low-metallicity IFMR of K12 with two IFMRs predicted by the detailed stellar evolution code MESA (Paxton et al. 2010) and with two IFMRs predicted by SeBa. The relation in age and mass in SeBa shows a strong upturn in age between WD masses with 0.54−0.53 M_⊙, after which the slopes become shallower again. This difference in slopes is due to two different evolution paths for low-mass main sequence stars. The higher mass stars follow a standard evolution path: they become WDs after the AGB phase. The lower mass stars on the other hand, lose their envelope on the RGB, whereafter they become WDs. In MESA this transition between these two evolution paths is implemented differently, yielding a more linear IFMR, which however is still steeper than the one inferred by K12. The observations seem to be consistent with all of these model lines. In fact, it will be challenging to observationally distinguish between the mass-age relations predicted by MESA and SeBa, since the difference between the two sets of lines is largest after a Hubble time. If on the other hand, ~0.51 M_⊙ single WDs will be found to follow the black solid line in Fig. 9, as some data seems to imply (see e.g. Table 1 of Renzini et al. 1996), there is a challenge for the theoretical modellers to explain how such low mass WDs can be formed by single stellar evolution within the age of the universe.

This comparison between MESA and SeBa was made using AMUSE (Portegies Zwart et al. 2009, 2013; Pelupessy et al. 2013). See also Renedo et al. (2010) for a comparison between theoretically and observationally determined IFMRs with different metallicities.

4. Conclusions

The easiest way to constrain the IMF with halo WDs is to determine the halo WD number density, because the normalization of the WDLF is linked one-to-one to the IMF. From a comparison between our derived halo WDLFs with the WDLF observed by RH11, we conclude that a Kroupa IMF (χ² = 2.29) is slightly preferred over a Salpeter IMF (χ² = 2.99). A top-heavy IMF (χ² = 5.74) clearly overpredicts the number of faint halo WDs. Due to large uncertainties on the normalization, it is not yet possible to completely rule out a non-standard IMF. However, also the shape of the WDLF corresponding to a top-heavy IMF has worse agreement with the WDLF observed by RH11 than those of the WDLFs corrsponding to models with a Kroupa or Salpeter IMF. Although most investigated halo models match the observed WDLF approximately equally well if we fix the normalization of our WDLF ($\mathrm{2.3}\lesssim {\mathit{\chi }}_{\min }^{\mathrm{2}}\lesssim \mathrm{2.7}$), none of the models comes close to a reduced χ² value of 1.

The exact number of halo WDs that Gaia can observe depends on how easily they can be distinguished from thin and thick disk WDs. In our standard model we find that Gaia will be able to detect approximately 1.5 × 10³ halo WDs, which is an order of magnitude more than the currently known number of halo WDs. Taking into account selection effects will probably not reduce this number by more than a factor of two. A wrong assumption on the mass function of unevolved stars has a stronger effect on the determined number density of halo WDs, but this will probably also be constrained by Gaia. If our assumptions are correct, the error bars on the observationally determined WDLF will become smaller with Gaia, at the part of the luminosity function that already has the smallest error bars (e.g. 5 ≲ M_bol ≲ 10), but also for the fainter luminosity bins. This implies that we might soon be able to start ruling out IMFs on the basis of their predicted WD number densities, especially in the M_bol ≳ 15 luminosity bins.

Since the effect of the SFH of halo stars is the strongest at the faint end of the WDLF, where we expect only a handful of WDs to be detected by Gaia, it will be observationally challenging to put strong constraints on this parameter in the near future. Although the differences are small, from the current observational constraints, we find that a model in which there was a burst of SF 13 Gyr ago (χ² = 2.29, ${\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{=}\mathrm{2.26}$) is slightly preferred over a burst of starformation 10 Gyr ago (χ² = 2.35, ${\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{=}\mathrm{2.43}$) and over continuous star formation 10 − 13 Gyr ago (χ² = 2.31, ${\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{=}\mathrm{2.37}$).

A determination of the masses of the brightest halo WDs can in theory be used to determine the age of the halo as suggested by K12. However, at this point only preliminary conclusions can be drawn since the observational uncertainties are large and the effect of binarity and/or metallicity can be underestimated. It would be useful to have more anchor points on the low-metallicity IFMR than the current single one.

With Gaia it will be possible to constrain for the first time the bright part of the field halo WDLF, where contributions from (unresolved) binary WDs are considerable (that this can be done in star clusters, and for singly evolved WDs was shown by Goldsbury et al. 2012). By determining the periods of WDs with masses below ≈0.5 M_⊙ (which can safely be assumed to be in binaries or to be the result of a binary merger) with Gaia follow-up observations, we can start to explore how binary stars with low metallicities evolve. In this paper, we do not vary any binary evolution parameters, but deviations from the predictions we made in the bottom panel of Fig. 6 are produced by different models of binary evolution.

It might be possible to put some constraints on the binary fraction by the number of pre-WDs that will be observed, although we expect this number to be small, since a large fraction of the pre-WD candidates will be indistinguishable from main sequence stars or giants. Furthermore, with the Bergeron models for WD cooling, pre-WDs are not expected to exist at all. However, pre-WDs can help us to constrain WD cooling models, because they are situated on an uncertain part of the WD cooling track, in the early phases of WD cooling.

If future observations on halo WDs go up to fainter magnitudes than Gaia can observe, we will be able to determine the validity of a top-heavy IMF or the Bergeron cooling models. In this respect, observations of the Large Synoptic Survey Telescope (LSST) will be very helpful (over 4 × 10⁵ halo WDs to r < 24.5; LSST Science Collaboration et al. 2009). To improve this study, WD cooling tracks and corresponding colours and magnitudes over the whole parameter range of WD masses and cooling times would be useful, for WDs with low-metallicity progenitors. In the near future, we will couple a semi-analytic model for galaxy formation with a binary population synthesis code and study how this affects the halo WD population.

Online material

Appendix A: Appendix A

Cartesian coordinates are related to spherical coordinates by $\begin{array}{ccc}& \mathit{x}& \mathrm{=}\mathit{r}\sin \mathit{\theta }\cos \mathit{\phi }\\ & \mathit{y}& \mathrm{=}\mathit{r}\sin \mathit{\theta }\sin \mathit{\phi }\\ & \mathit{z}& \mathrm{=}\mathit{r}\cos \mathit{\theta }\end{array}$with radius r, polar angle θ and azimuth angle φ. The Sun is assumed to be at position (x,y,z)_⊙ = (r₀,0,0), or equivalently at (r,θ,φ)_⊙ = (r₀,π/ 2,0). We define the primed coordinates $\begin{array}{ccc}& {\mathit{x}}^{\mathrm{\prime }}& \mathrm{=}{\mathit{r}}^{\mathrm{\prime }}\sin {\mathit{\theta }}^{\mathrm{\prime }}\cos \mathit{\phi }\mathrm{\equiv }\mathit{x}\\ & {\mathit{y}}^{\mathrm{\prime }}& \mathrm{=}{\mathit{r}}^{\mathrm{\prime }}\sin {\mathit{\theta }}^{\mathrm{\prime }}\sin \mathit{\phi }\mathrm{\equiv }\mathit{y}\\ & {\mathit{z}}^{\mathrm{\prime }}& \mathrm{=}{\mathit{r}}^{\mathrm{\prime }}\cos {\mathit{\theta }}^{\mathrm{\prime }}\mathrm{\equiv }\mathit{q}\mathit{z}\end{array}$such that the local halo density (Eq. (5)) can be expressed independent of a polar angle and azimuth angle: $\mathit{\rho }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{y}}^{\mathrm{\prime }}\mathit{,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{\equiv }\mathit{\rho }\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\frac{{\mathit{r}}^{\mathrm{\prime }}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathit{n}}\mathrm{·}$(A.7)We note that in the Galactic plane, z = 0, thus primed radius r′ = r and the primed polar angle θ′ = θ = ^π/ 2. At the Galactic pole, θ′ = θ = 0, and r′ = z′ = qz = qr. In all other cases, the relation between the r′, θ′ and their spherical equivalents is given by $\begin{array}{ccc}& {\mathit{\theta }}^{\mathrm{\prime }}& \\ & {\mathit{r}}^{\mathrm{\prime }}& \end{array}$Since we assume an oblate stellar halo (q < 1), it follows from Eqs. (A.8) and (A.9) that θ′ ≥ θ and r′ ≤ r for any given point in the spheroid. Because we want a sphere with radius ξ around the Sun to be contained in our simulated area, we set the boundary conditions, $\begin{array}{ccc}{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\xi }& \mathrm{\le }\mathit{r}& \mathrm{\le }{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\xi }\\ \mathit{\delta }& \mathrm{\le }\mathit{\theta }& \mathrm{\le }\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{-}\mathit{\delta }\\ \mathrm{-}\mathit{ϵ}& \mathrm{\le }\mathit{\phi }& \mathrm{\le }\mathit{ϵ,}\end{array}$with δ ≤ arctan(^r₀q/ξ) and ϵ ≤ arctan(^ξ/r₀). These set the limits of integration in our determination of the stellar halo mass: $\mathit{M}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{0}}^{\mathit{n}}}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\xi }}^{{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\xi }}{\mathit{r}}^{\mathit{n}\mathrm{+}\mathrm{2}}\mathrm{d}\mathit{r}{\mathrm{\int }}_{\mathit{\delta }}^{\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{-}\mathit{\delta }}{\left(\frac{{\cos }^{\mathrm{2}}\mathit{\theta }}{{\mathit{q}}^{\mathrm{2}}}\mathrm{+}{\sin }^{\mathrm{2}}\mathit{\theta }\right)}^{\mathit{n}\mathit{/}\mathrm{2}}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }{\mathrm{\int }}_{\mathrm{-}\mathit{ϵ}}^{\mathit{ϵ}}\mathrm{d}\mathit{\phi }\mathit{.}$(A.13)In order to solve the integral over θ, we now first make an estimation of δ. With the assumed values of ξ, q and r₀ mentioned in the main text, we find δ ≤ 0.334 π. Thus, we take δ = π/ 3. The integral over θ can now be expressed as the hypergeometric function ${\mathrm{2}}_{}{\mathit{F}}_{\mathrm{1}}{}^{\mathrm{(}}\mathrm{1}^{\mathrm{1}}{}^{\mathrm{1}}\mathit{/}\mathrm{2}\mathit{,}\mathrm{-}{}^{\mathit{n}}\mathit{/}\mathrm{2}\mathrm{;}{}^{\mathrm{3}}\mathit{/}\mathrm{2}\mathrm{;}{}^{\mathrm{1}}\mathit{/}\mathrm{4}\mathrm{-}{}^{\mathrm{1}}\mathit{/}{\mathrm{4}{\mathit{q}}^{\mathrm{2}}}^{\mathrm{)}}$. Again with q = 0.64 and n = −2.8 for consistency with Jurić et al. (2008), we find ${\mathrm{2}}_{}{\mathit{F}}_{\mathrm{1}}\left(\mathrm{0.5}\mathit{,}\mathrm{1.4}\mathrm{;}\mathrm{1.5}\mathrm{;}\mathrm{-}\mathrm{0.36}\right)\mathrm{=}\mathrm{0.866}$. Because this value of n ≠ −3, the integral over r can also be evaluated: $\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{0}}^{\mathit{n}}}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\xi }}^{{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\xi }}{\mathit{r}}^{\mathit{n}\mathrm{+}\mathrm{2}}\mathrm{d}\mathit{r}\mathrm{=}\frac{\mathrm{(}{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\xi }{\mathrm{)}}^{\mathit{n}\mathrm{+}\mathrm{3}}\mathrm{-}\mathrm{(}{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\xi }{\mathrm{)}}^{\mathit{n}\mathrm{+}\mathrm{3}}}{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{3}\mathrm{)}{\mathit{r}}_{\mathrm{0}}^{\mathit{n}}}\mathrm{=}\mathrm{3.91}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathrm{pc}}^{\mathrm{3}}\mathit{.}$(A.14)The integral over φ yields 2ϵ, thus after choosing ϵ = arctan(^ξ/r₀) this reads 2arctan(^ξ/r₀) = 0.707. The multiplication of an assumed value of ρ₀ = 1.5 × 10^-4 M_⊙ pc^-3 (Fuchs & Jahreiß 1998) with these three integrals gives M_unev = 3.6 × 10⁷ M_⊙.

Appendix B: Appendix B

In case φ(m) is a single power law function between the upper and lower mass boundary of unevolved stars in our simulation box m_high,unev and m_low,unev, the total mass in unevolved stars $\frac{{\mathit{M}}_{\mathrm{unev}}}{{\mathit{M}}_{\mathrm{\odot }}}\mathrm{=}{\mathrm{\int }}_{{\mathit{m}}_{\mathrm{low}\mathit{,}\mathrm{unev}}}^{{\mathit{m}}_{\mathrm{high}\mathit{,}\mathrm{unev}}}\mathit{A}{\mathit{m}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{unev}}}\mathrm{d}\mathit{m}\mathrm{=}\frac{\mathit{A}\left({\mathit{m}}_{\mathrm{high}\mathit{,}\mathrm{unev}}^{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{unev}}}\mathrm{-}{\mathit{m}}_{\mathrm{low}\mathit{,}\mathrm{unev}}^{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{unev}}}\right)}{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{unev}}}\mathrm{·}$(B.1)Given the mass in unevoloved stars M_unev which was derived in Appendix A, γ_unev = −1, m_high,unev = 0.8 and m_low,unev = 0.1, this results in a normalization constant belonging to the lower limit on the number of unevolved (single) stars N_unev in our simulation box A_lower = 1.1 × 10⁸. When substituted into Eq. (4), this yields ${\mathit{N}}_{\mathrm{unev}}\mathit{>}{\mathit{A}}_{\mathrm{lower}}\mathrm{(}{\mathit{m}}_{\mathrm{high}\mathit{,}\mathrm{unev}}\mathrm{-}{\mathit{m}}_{\mathrm{low}\mathit{,}\mathrm{unev}}\mathrm{)}\mathrm{=}\mathrm{8.0}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\mathit{.}$(B.2)We derive an upper limit on the number of evolved stars N_ev in our simulation box, for the three different IMFs that we investigate in this paper by determining their normalization constants from the IMF at m_high,unev. For example, writing the normalization constant for the upper limit on the number of evolved stars in case of a Kroupa IMF as B_upper, the relation φ(m_high,unev) = A_lower = B_upper (m_high,unev)^-2.2 leads to B_upper = 7.0 × 10⁷, from which follows ${\mathit{N}}_{\mathrm{ev}\mathit{,}\mathrm{Kroupa}}\mathit{<}{\mathit{B}}_{\mathrm{upper}}\mathrm{·}{\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Kroupa}}\mathrm{=}\mathrm{5.9}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\mathit{,}$(B.3)where ${\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Kroupa}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0.8}}^{\mathrm{1.0}}{\mathit{m}}^{-2.2}\mathrm{d}\mathit{m}\mathrm{+}{\mathrm{\int }}_{\mathrm{1.0}}^{\mathrm{100}}{\mathit{m}}^{-2.7}\mathrm{d}\mathit{m}\mathrm{=}\mathrm{0.84.}$(B.4)To obtain actual numbers instead of an upper limit, we assume that the low-mass part of the IMF is correctly given by Eq. (2), with normalization constant B, $\frac{{\mathit{M}}_{\mathrm{unev}}}{{\mathit{M}}_{\mathrm{\odot }}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0.1}}^{\mathrm{0.5}}\frac{\mathrm{35}}{\mathrm{19}}\mathit{B}{\mathit{m}}^{-0.3}\mathrm{d}\mathit{m}\mathrm{+}{\mathrm{\int }}_{\mathrm{0.5}}^{\mathrm{0.8}}\mathit{B}{\mathit{m}}^{-1.2}\mathrm{d}\mathit{m}\mathrm{=}\mathrm{1.6}\mathit{B,}$(B.5)again using the calculated total mass in unevolved stars M_unev = 3.6 × 10⁷ M_⊙, we find B = 2.2 × 10⁷. Now because ${\mathit{I}}_{\mathrm{unev}\mathit{,}\mathrm{Kroupa}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0.1}}^{\mathrm{0.5}}\frac{\mathrm{35}}{\mathrm{19}}{\mathit{m}}^{-1.3}\mathrm{d}\mathit{m}\mathrm{+}{\mathrm{\int }}_{\mathrm{0.5}}^{\mathrm{0.8}}{\mathit{m}}^{-2.2}\mathrm{d}\mathit{m}\mathrm{=}\mathrm{5.5}$(B.6)we find $\begin{array}{ccc}& & {\mathit{N}}_{\mathrm{unev}\mathit{,}\mathrm{Kroupa}}\mathrm{=}\mathit{B}\mathrm{·}{\mathit{I}}_{\mathrm{unev}\mathit{,}\mathrm{Kroupa}}\mathrm{=}\mathrm{1.2}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}\\ & & {\mathit{N}}_{\mathrm{ev}\mathit{,}\mathrm{Kroupa}}\mathrm{=}\mathit{B}\mathrm{·}{\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Kroupa}}\mathrm{=}\mathrm{1.9}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\mathit{.}\end{array}$Assuming that the Salpeter IMF holds for masses m> 0.8 results in the same way into an upper limit on the number of evolved stars, whereas assuming that it is for the entire mass range 0.1 <m < 100 gives the expected number of evolved stars. Since ${\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0.8}}^{\mathrm{100}}{\mathit{m}}^{-2.35}\mathrm{d}\mathit{m}\mathrm{=}\mathrm{1.00}\mathit{,}$(B.9)the upper limit on the number of evolved stars in the case of a Salpeter IMF immediately follows from the normalization constant ${\mathit{C}}_{\mathrm{upper}}\mathrm{=}{\mathit{A}}_{\mathrm{lower}}\mathit{/}{\mathit{m}}_{\mathrm{high}\mathit{,}\mathrm{unev}}^{-2.35}\mathrm{=}\mathrm{6.7}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}$, ${\mathit{N}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}}\mathit{<}{\mathit{C}}_{\mathrm{upper}}\mathrm{·}{\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}\mathrm{6.7}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\mathit{.}$(B.10)The expected number of stars in our simulation box if the low-mass part of the mass function is also Salpeter $\begin{array}{ccc}& & {\mathit{N}}_{\mathrm{unev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}\mathit{C}\mathrm{·}{\mathit{I}}_{\mathrm{unev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}\mathrm{1.7}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}\\ & & {\mathit{N}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}\mathit{C}\mathrm{·}{\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}\mathrm{1.1}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\end{array}$with $\frac{{\mathit{M}}_{\mathrm{unev}}}{{\mathit{M}}_{\mathrm{\odot }}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0.1}}^{\mathrm{0.8}}\mathit{C}{\mathit{m}}^{-1.35}\mathrm{d}\mathit{m}\mathrm{=}\mathrm{3.3}\mathit{C,}$(B.13)thus C = 1.1 × 10⁷, and ${\mathit{I}}_{\mathrm{unev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0.1}}^{\mathrm{0.8}}{\mathit{m}}^{-2.35}\mathrm{d}\mathit{m}\mathrm{=}\mathrm{15.6.}$(B.14)Finally, for the top-heay IMF we derive the normalization constants for the Komiya IMF (indicated by the letter D) and the Salpeter IMF (indicated by the letter E) simultaneously, using the MDF of the halo described by An et al. (2013), who studied halo main-sequence stars with masses between 0.65  M_⊙ and 0.75  M_⊙ in the Sloan Digital Sky Survey. These authors found that the halo can be described by a two-component model, with 24% of the stars belonging to a low-metallicity population with a peak at [Fe/H] = − 2.33 (i.e. their calibration model). If this population of low-metallicity stars is born according to a Komiya IMF, we have ${\mathrm{\int }}_{\mathrm{0.65}}^{\mathrm{0.75}}\mathit{D}\exp \left[\mathrm{-}\frac{{\log }_{\mathrm{10}}^{\mathrm{2}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\frac{\mathrm{d}\mathit{m}}{\mathit{m}}\mathrm{=}\frac{\mathrm{0.24}}{\mathrm{0.76}}{\mathrm{\int }}_{\mathrm{0.65}}^{\mathrm{0.75}}\mathit{E}{\mathit{m}}^{-2.35}\mathrm{d}\mathit{m,}$(B.15)which holds for and D and E, as well as for D_upper and E_upper. The normalization constants for the upper limit on the number of evolved stars in case of a top-heavy IMF follow again from $\begin{array}{ccc}& \mathit{\phi }\mathrm{\left(}{\mathit{m}}_{\mathrm{high}\mathit{,}\mathrm{unev}}\mathrm{\right)}& \mathrm{=}\frac{{\mathit{D}}_{\mathrm{upper}}}{{\mathit{m}}_{\mathrm{high}\mathit{,}\mathrm{unev}}}\exp \left[\mathrm{-}\frac{{\log }_{\mathrm{10}}^{\mathrm{2}}\mathrm{(}{\mathit{m}}_{\mathrm{high}\mathit{,}\mathrm{unev}}\mathit{/}\mathit{\mu }\mathrm{)}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\\ & & \mathrm{+}{\mathit{E}}_{\mathrm{upper}}\mathrm{(}{\mathit{m}}_{\mathrm{high}\mathit{,}\mathrm{unev}}{\mathrm{)}}^{-2.35}\mathrm{=}{\mathit{A}}_{\mathrm{lower}}\mathit{.}\end{array}$(B.16)

From the standard integral $\mathrm{\int }\exp \left[\mathrm{-}\frac{{\log }_{\mathrm{10}}^{\mathrm{2}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\frac{\mathrm{d}\mathit{m}}{\mathit{m}}\mathrm{=}\frac{\sqrt{\mathit{\pi }\mathit{/}\mathrm{2}}\mathit{\sigma }\mathrm{erf}\left[{\log }_{\mathrm{10}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}\mathit{/}\sqrt{\mathrm{2}}\mathit{\sigma }\right]}{{\log }_{\mathrm{10}}\mathit{e}}$(B.17)it now follows that D_upper = 1.4 × 10⁹ and E_upper = 4.3 × 10⁷. Consequently, the number of evolved stars $\begin{array}{ccc}& & {\mathit{N}}_{\mathrm{ev}\mathit{,}\mathrm{Komiya}}\mathit{<}{\mathit{D}}_{\mathrm{upper}}\mathrm{·}{\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Komiya}}\mathrm{=}\mathrm{3.3}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}\\ & & {\mathit{N}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}\mathrm{(}\mathrm{top}\mathrm{-}\mathrm{heavy}\mathrm{)}}\mathit{<}{\mathit{E}}_{\mathrm{upper}}\mathrm{·}{\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}\mathrm{4.3}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\end{array}$with ${\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Komiya}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0.8}}^{\mathrm{100}}\exp \left[\mathrm{-}\frac{{\log }_{\mathrm{10}}^{\mathrm{2}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\frac{\mathrm{d}\mathit{m}}{\mathit{m}}\mathrm{=}\mathrm{2.29.}$(B.20)If the suggested top-heavy IMF holds in the low-mass regime, $\begin{array}{ccc}\frac{{\mathit{M}}_{\mathrm{unev}}}{{\mathit{M}}_{\mathrm{\odot }}}& & \mathrm{=}{\mathrm{\int }}_{\mathrm{0.1}}^{\mathrm{0.8}}\mathit{D}\exp \left[\mathrm{-}\frac{{\log }_{\mathrm{10}}^{\mathrm{2}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\mathrm{d}\mathit{m}\mathrm{+}{\mathrm{\int }}_{\mathrm{0.1}}^{\mathrm{0.8}}\mathit{E}{\mathit{m}}^{-1.35}\mathrm{d}\mathit{m}\\ & & \end{array}$(B.21)where we used the standard integral: $\begin{array}{ccc}\mathrm{\int }\exp \left[\mathrm{-}\frac{{\log }_{\mathrm{10}}^{\mathrm{2}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\mathrm{d}\mathit{m}& \mathrm{=}& \\ & & \mathrm{-}\frac{\sqrt{\mathit{\pi }\mathit{/}\mathrm{2}}\mathit{\mu \sigma }\exp \left(\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}{\log }_{\mathrm{10}}^{\mathrm{2}}\mathit{e}}\right)}{{\log }_{\mathrm{10}}\mathit{e}}\mathrm{erf}\left[\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}{\log }_{\mathrm{10}}\mathit{e}{\log }_{\mathrm{10}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}}{\sqrt{\mathrm{2}}\mathit{\sigma }{\log }_{\mathrm{10}}\mathit{e}}\right]\mathit{.}\end{array}$(B.22)Combining Eqs. (B.15) and (B.21), we find D = 3.4 × 10⁸ and E = 1.0 × 10⁷, as well as $\begin{array}{ccc}& & {\mathit{N}}_{\mathrm{unev}\mathit{,}\mathrm{Komiya}}\mathrm{=}\mathit{D}\mathrm{·}{\mathit{I}}_{\mathrm{unev}\mathit{,}\mathrm{Komiya}}\mathrm{=}\mathrm{2.4}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}\\ & & {\mathit{N}}_{\mathrm{unev}\mathit{,}\mathrm{Salpeter}\mathrm{(}\mathrm{top}\mathrm{-}\mathrm{heavy}\mathrm{)}}\mathrm{=}\mathit{E}\mathrm{·}{\mathit{I}}_{\mathrm{unev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}\mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}\\ & & {\mathit{N}}_{\mathrm{ev}\mathit{,}\mathrm{Komiya}}\mathrm{=}\mathit{D}\mathrm{·}{\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Komiya}}\mathrm{=}\mathrm{7.9}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}\\ & & {\mathit{N}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}\mathrm{(}\mathrm{top}\mathrm{-}\mathrm{heavy}\mathrm{)}}\mathrm{=}\mathit{E}\mathrm{·}{\mathit{I}}_{\mathrm{ev}\mathit{,}\mathrm{Salpeter}}\mathrm{=}\mathrm{1.0}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\mathit{,}\end{array}$where ${\mathit{I}}_{\mathrm{unev}\mathit{,}\mathrm{Komiya}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0.1}}^{\mathrm{0.8}}\exp \left[\mathrm{-}\frac{{\log }_{\mathrm{10}}^{\mathrm{2}}\mathrm{(}\mathit{m}\mathit{/}\mathit{\mu }\mathrm{)}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\frac{\mathrm{d}\mathit{m}}{\mathit{m}}\mathrm{=}\mathrm{7.0}\mathrm{\times }{\mathrm{10}}^{-3}\mathit{.}$(B.27)

Acknowledgments

We thank Else Starkenburg as well as the anonymous referee for valuable comments. We thank the Netherlands Research School for Astronomy (NOVA) and the Netherlands Organisation for Scientific Research (NWO) for financial support (#639.073.803 [VICI] and #614.061.608 [AMUSE]). A.H. acknowledges financial support from ERC StG 240271, Galactica.

References

All Tables

All Figures

	Fig. 1 Simulation box containing a sphere with radius ξ = 2.95 kpc centred at the position of the Sun (x,y,z)⊙ = (8.0,0,0) kpc, n = − 2.8 and q = 0.64 (left panel) and a density map of the simulation box (right panel).
In the text

Fig. 2 Comparison of cooling tracks of 0.2  M⊙ (dashed lines), 0.5  M⊙ (solid lines) and 1.2  M⊙ (dotted lines) WDs, between the models from the Althaus group (dark blue lines) and the Bergeron group (light blue lines). The 3 lines that represent cooling tracks from the Althaus models correspond to 3 WDs with different core compositions: He for the 0.2  M⊙ WD, CO for the 0.5  M⊙ WD and ONe for the 1.2  M⊙ WD, whereas the 3 lines from the Bergeron models all correspond to WDs with a CO core. The age of the universe is indicated with a vertical thin black line.

In the text

Fig. 3 From left to right: constructed grids of He, CO and ONe WD cooling in the Althaus models. The width of the panels corresponds to the maximum possible mass range per panel in our simulation, see Table 2. Vertical (horizontal) dashed lines indicate masses (cooling times) beyond which luminosities are obtained from extrapolation of the used cooling tracks. Dotted lines indicate the boundaries beyond which MV and MI magnitudes are obtained by extrapolation. For ONe WDs these regions completely overlap, for CO WDs the lower mass boundary partly overlaps, which is indicated by the light-grey-dark-grey dashed line. Scattered points indicate halo WDs that are observable with Gaia in our standard model. An explanation of the symbols is given Sect. 3.3.

In the text

Fig. 4 Build-up of the WDLF in our standard halo model (50% binaries). All WDs are included in the black solid line (the total WDLF). The lower solid line shows the contribution from unresolved binary WDs. The dashed, dot-dashed and dotted lines show the contributions from CO, He and ONe WDs respectively. The light blue line with error bars is the halo WDLF determined from the SuperCOSMOS Sky Survey (Rowell & Hambly 2011), shown for comparison.

In the text

	Fig. 5 Comparison of halo WD luminosity functions corresponding to different assumptions about the IMF (top left), WD cooling (top right), binary fraction (bottom left) and SFH of the halo (bottom right). The WDLF corresponding to our standard model, indicated by the blue solid line in all panels, is constructed using the 1.2 × 107 WDs in our simulation box (see Table 1).
In the text

Fig. 6 Properties of unresolved binary WD pairs with G < 20 for our standard model. Shown are double He WDs (circles), double CO WDs (squares), CO+He systems in which the He WD is the brightest (upward pointing triangles), He+CO systems in which the CO WD is the brightest (downward pointing traingles), pre-WDs (filled stars), and “other” stars, which are pre-WDs that are indistinguishable from MS stars or giants. After the label discriptions in the legend, the total number of WD binaries of that particular kind is given. Due to the low number of halo WDs with G < 20 we expect to find in the Milky Way, there is some statistical noise in this figure.

In the text

Fig. 7 Properties of all 1.3 × 106 unresolved binary WD pairs in our simulation box (maximum magnitude G = 35), for our standard model. Contour lines mark the regions in which 33%, 67% and 95% of the binaries are located. Higher density regions are given darker colours. The 1% binaries in the tail of the distribution are indicated by the scattered points. Compared to Fig. 6, this figure hardly has statistical noise because of the large number of halo WDs we expect to find in the Milky Way halo if we would go up to magnitude G = 35.

In the text

Fig. 8 Distances and bolometric magnitudes of the 1.5 × 103 halo WDs that can be observed with Gaia, for our standard model. Top panel: distribution of their distances. Right panel: distribution of their bolometric magnitudes, which gives an idea of the statistical errors that are to be expected per luminosity bin of the WDLF for Gaia. The yellow stars, labelled “other”, are pre-WDs that are indistinguishable from MS stars or giants. These are not included in the projected distribution histograms.

In the text

Fig. 9 Age of the brightest single halo WDs (and thus of the Galactic halo) as a function of WD mass. Since for single stars every age corresponds one-to-one to an initial stellar mass and WD mass, this is an alternative representation of the IFMR. The IFMRs for single stars predicted by SeBa (light blue curves) and the those predicted by MESA (dark blue curves) both differ from the black solid straight line through the data points with error bars of Kalirai (2012). The dashed lines correspond to the metallicity value that we used as our standard for the Milky Way halo in this study (Z = 0.001), while the dotted lines correspond to half that metallicity value (Z = 0.0005).

In the text