# Local stellar formation history from the 40 pc white dwarf sample

E. Cukanovaite1 111E-mail<EMAIL_ADDRESS>P.-E. Tremblay1, S.
Toonen2, K.D. Temmink3, Christopher J. Manser4, M.W. O’Brien1 and J. McCleery1
1 Department of Physics, University of Warwick, Coventry CV4 7AL, UK
2 Anton Pannekoek Institute for Astronomy, University of Amsterdam, 1090 GE
Amsterdam, The Netherlands
3 Department of Astrophysics/IMAPP, Radboud University Nijmegen, P.O. Box
9010, 6500 GL Nijmegen, The Netherlands
4 Astrophysics Group, Department of Physics, Imperial College London, Prince
Consort Rd, London, SW7 2AZ, UK

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We derive the local stellar formation history from the Gaia-defined 40 pc
white dwarf sample. This is currently the largest volume-complete sample of
white dwarfs for which spectroscopy is available, allowing for classification
of the chemical abundances at the photosphere, and subsequently accurate
determination of the atmospheric parameters. We create a population synthesis
model and show that a uniform stellar formation history for the last $\approx
10.5$ Gyr provides a satisfactory fit to the observed distribution of absolute
Gaia $G$ magnitudes. To test the robustness of our derivation, we vary various
assumptions in the population synthesis model, including the initial mass
function, initial-to-final mass relation, kinematic evolution, binary fraction
and white dwarf cooling timescales. From these tests, we conclude that the
assumptions in our model have an insignificant effect on the derived relative
stellar formation rate as a function of look-back time. However, the onset of
stellar formation (age of Galactic disc) is sensitive to a variety of input
parameters including the white dwarf cooling models. Our derived stellar
formation history gives a much better fit to the absolute Gaia $G$ magnitudes
than most previous studies.

###### keywords:

white dwarfs – stars: luminosity function – stars: evolution

††pubyear: 2022††pagerange: Local stellar formation history from the 40 pc
white dwarf sample–Local stellar formation history from the 40 pc white dwarf
sample

##  1 Introduction

Clues of the past and future of the universe can be found by studying galaxies
at different redshifts and therefore at different epochs of time (see, for
example, Speagle et al. 2014; Somerville & Davé 2015; Frebel & Norris 2015;
Kruijssen et al. 2019 and references therein). Our own Milky Way is a perfect
representative of spiral galaxies in the low-redshift universe (see, for
example, Helmi 2008; Hou & Han 2014; Helmi 2020 and references therein)
offering the additional advantage of being the only galaxy that can be studied
from within. By examining the evolution of the individual stars in the Milky
Way, we can learn how the baryonic matter of galaxies has developed over time.

The overall structure of the Milky Way is thought to be made up of four main
baryonic components: the bulge and the bar, the thin disc, the thick disc, and
the halo (Helmi, 2020). It is agreed that the Milky Way follows the
$\Lambda$CDM cosmological model (White & Rees, 1978; Davis et al., 1985; De
Lucia et al., 2006; Planck Collaboration et al., 2016; Bullock & Boylan-
Kolchin, 2017), where galaxies form inside dark matter halos, attracting
baryonic matter over time. The baryonic matter then cools and collapses
towards the centre of the dark matter halo, forming a disc due to conservation
of angular momentum (Mo et al., 1998; Cole et al., 2000; Kereš et al., 2005).
Most importantly, the galaxy also grows via mergers with neighbouring galaxies
(Odenkirchen et al., 2001; Ibata et al., 2002; Niederste-Ostholt et al., 2010;
Besla et al., 2010; Belokurov, 2013; Grillmair & Carlin, 2016; Helmi, 2020),
with the first point of contact being the outer halo of the Milky Way.
Nevertheless, the mergers have an impact on all the components of the Galaxy.
These mergers leave imprints in the form of Galactic streams, which have been
observed and studied extensively (Belokurov, 2013; Grillmair & Carlin, 2016;
Helmi, 2020).

It is accepted that all the components of the Galaxy are interlinked in their
evolution, due to a single Galactic potential (Guedes et al., 2013) and radial
migration of stars as they get older (Sellwood & Binney, 2002; Siebert et al.,
2008; Bond et al., 2010; Minchev & Famaey, 2010; Minchev et al., 2015). Thus,
some studies have trouble separating the different components (Norris & Ryan,
1991; Bovy et al., 2012a; Bovy et al., 2012b; Hayden et al., 2017). It is now
accepted that kinematically the thin and thick disc stars are hard to
differentiate, yet, the two components can be separated by stellar metallicity
(Minchev et al., 2015; Kawata & Chiappini, 2016; Helmi, 2020), with the
members of the thick disc being more metal-poor and therefore older (Gilmore
et al., 1989; Haywood et al., 2013; Hayden et al., 2015). Similar issues have
recently been uncovered in trying to separate the thick disc from the inner
halo, as both are thought to have been significantly affected by the Gaia-
Enceladus merger, and are observed to overlap in the Solar neighbourhood
(Koppelman et al., 2018; Belokurov et al., 2018; Haywood et al., 2018; Helmi
et al., 2018; Gallart et al., 2019; Myeong et al., 2019; Helmi, 2020).
However, it can be summarised that the oldest stars are found in the halo,
which was formed around 10-12 Gyr ago (Salaris & Weiss, 2002; Jofré & Weiss,
2011; Kilic et al., 2017; Gallart et al., 2019). On the other hand, the thin
disc is the youngest, formed around 8 Gyr ago (Jimenez et al., 1998; Liu &
Chaboyer, 2000; del Peloso et al., 2005a, b; Haywood et al., 2013; Kilic et
al., 2017; Tononi et al., 2019). Ultimately, it is apparent that the evolution
of the Milky Way is complex and is not yet fully understood, but has recently
been and continues to be aided by the data coming from the Gaia satellite
(Gaia Collaboration et al., 2016a, b, 2018, 2021).

The evolution of the Milky Way can be summarised in a stellar formation
history, which denotes how many stars were formed during a particular period.
Traditionally, its derivation has relied on methods involving main sequence
stars, but other types of stars can also be used. Some example methods include
nucleocosmochronometry, isochrone fitting, and asteroseismology (see, for
example, the review of Soderblom 2010 and references therein). An independent
and alternative method is based on white dwarf stars (see, for example, Winget
et al. 1987; Rowell 2013; Tremblay et al. 2014; Kilic et al. 2017; Isern 2019;
Fantin et al. 2019). These stars do not undergo any significant nuclear fusion
in their cores and instead cool with age via the gravothermal process,
providing accurate age estimates (Althaus et al., 2010a; Rebassa-Mansergas et
al., 2015; Heintz et al., 2022). Additionally, white dwarfs encompass
information over the entirety of the Milky Way’s history, with the oldest
known white dwarfs having cooling ages of around 10 Gyr, meaning that they
have come from main sequence progenitors that are even older (Fontaine et al.,
2001).

The first established method of using white dwarfs to determine the stellar
formation history involved their luminosity function (Winget et al., 1987),
where a sharp drop-off is observed at dim magnitudes. This was understood as
the finite age of the Milky Way and it being smaller than the time needed for
a white dwarf to cool down to magnitudes dimmer than the drop-off in the
luminosity function. This method has been applied both to the Galaxy as a
whole (Winget et al., 1987; Garcia-Berro et al., 1988; Bessell & Stringfellow,
1993) and to the individual Galactic components (see, for example, Mochkovitch
et al. 1990; Isern et al. 1998; García-Berro et al. 2010; Kilic et al. 2017;
Torres et al. 2019; Fantin et al. 2019, 2021; Torres et al. 2021 and reference
therein). Furthermore, these and other studies have also shown that the bright
end of the white dwarf luminosity function is also highly sensitive to the
stellar formation history (Noh & Scalo, 1990; García-Berro & Oswalt, 2016).
Thus, the white dwarf luminosity function tracks the entirety of the Galactic
stellar formation history and can be manipulated to get more accurate results.
For example, Isern (2019) and references therein showed that massive white
dwarfs can be used to find stellar formation histories without needing to
invoke uncertainties from the main sequence lifetimes.

For either white dwarf or main sequence methods, the derivation of an accurate
stellar formation history is heavily impeded by various observational biases.
In particular, to overcome incompleteness of observational samples, population
synthesis codes have to include complex bias corrections, which, for example,
can be functions of survey pointing, magnitude and colour. With the advent of
Gaia and its data releases, the completeness of the night sky has increased
significantly (Gaia Collaboration et al., 2016a, b, 2018, 2021) with
unprecedented precision in astrometry (Lindegren et al., 2021b, a). One such
tremendous improvement involved the careful cataloguing of white dwarfs in
Gaia Data Release 2 (DR2, Gentile Fusillo et al. 2019) and early Data Release
3 (eDR3, Gentile Fusillo et al. 2021).

In this paper we aim to derive the local stellar formation history using a
subset of 1083 Gaia white dwarfs within 40 pc of the Sun. This subset is
estimated to have a Gaia completeness of around 96% (Tremblay et al., 2020;
McCleery et al., 2020). Furthermore, McCleery et al. (2020), Tremblay et al.
(2020) and O’Brien et al. (2023) have compiled the atmospheric compositions
for nearly all white dwarfs ($>$96%) in that sub-sample. They used Gaia-
independent spectroscopic observations, making the 40 pc sample the largest,
most complete volume of spectroscopically confirmed white dwarfs to date. This
is particularly important as the atmospheric composition of a white dwarf
directly affects its temperature, radius and mass derived from the Gaia data
(Bergeron et al., 2019), leading to errors of up to 20% in these parameters if
the atmospheric composition is instead assumed. Another advantage of this
sample compared to the larger volume samples is the negligible amount of
interstellar dust, allowing for a direct comparison of observed and predicted
Gaia $G$ magnitudes. Additionally, for larger samples the Gaia errors become
larger, making the parameters provided by Gaia less accurate, as well as
making the sampled volume harder to define.

In Sect. 2 we introduce the 40 pc Gaia sample of white dwarfs. Sect. 3
describes the simulation constructed to model the local white dwarf sample.
Sect. 4 describes the method used to derive the local stellar formation
history. We compare our results to other studies in Sect. 5. In Sect. 6 we
discuss the various uncertainties associated with assumptions made in the
simulation. We discuss the results and conclude in Sect. 7.

##  2 The 40 pc sample of white dwarfs

To derive the local stellar formation history, we use the northern 40 pc
sample of spectroscopically confirmed white dwarfs curated by McCleery et al.
(2020) and the southern 40 pc sample curated by O’Brien et al. (2023). The
northern catalogue was originally published based on DR2, but has recently
been updated with Gaia eDR3 (Gentile Fusillo et al.
2021)222https://cygnus.astro.warwick.ac.uk/phrgwr/40pcTables/index.html. In
total, the 40 pc sample contains 1083 white dwarfs: 626 have hydrogen-
dominated atmospheres; 266 have helium-dominated atmospheres; 9 are likely
helium-rich DC white dwarfs with strong CIA absorption. The remaining 182
white dwarfs have spectra that show no significant absorption lines and have
Gaia photometric $T_{\rm{eff}}$ below $\approx$ 5000 K, where both He and H
lines become invisible. Therefore, these cool white dwarfs have unconstrained
atmospheric compositions (see McCleery et al. 2020 for more detail).

Given the precision of the Gaia data, the observed absolute $G$ magnitude
provides a robust independent parameter for the determination of the local
stellar formation history. The absolute magnitude is calculated from the
distance modulus equation, where the distance is assumed to be equal to the
inverse of parallax. This is not a correct assumption for parallaxes with
large errors because of the data processing done by Gaia (Bailer-Jones et al.,
2018, 2021). However, as our sample is local and has small parallax errors,
this assumption is valid, and as such we utilise it. The distance modulus
equation also relies on the observed apparent $G$ magnitude, which given the
locality of our sample, we do not deredden, as dust becomes only important for
distances larger than around 50 pc (Gentile Fusillo et al., 2021).

Our approach is similar to earlier studies deriving the local stellar
formation history from the white dwarf luminosity function, but does not rely
on the determination of individual atmospheric parameters, the effective
temperature ($T_{\rm eff}$) and surface gravity, which were necessary in the
pre-Gaia era to infer the distance, and hence luminosity. In particular, our
method only depends on the absolute optical Gaia $G$ flux and does not use
Gaia colours (or Gaia inferred $T_{\rm eff}$), which are much more sensitive
to the details of the modelled atmosphere, and thus the opacities.
Nevertheless, any significant offset between observed and predicted Gaia
colours, such as that seen for white dwarfs with $T_{\rm eff}\lesssim 5000$ K
which suffer from a low-mass problem (Hollands et al., 2018; McCleery et al.,
2020), could have a smaller repercussion on the absolute Gaia $G$ flux.

In Fig. 1 the mass-$T_{\rm eff}$ distribution of the 40 pc sample is shown as
black dots. These parameters are derived using the method described in
McCleery et al. (2020), which is based on eDR3 colours and the appropriate
atmospheric composition given the identified spectral type.

There is a lower limit to the white dwarf mass for single star evolution,
since the creation of a white dwarf with a mass smaller than around 0.54 M⊙
would take longer than the current age of the Universe (see, for example,
Cummings et al. 2016). Thus, the very low observed masses in Fig. 1 are caused
by binary systems. In the case of unresolved double white dwarf systems, the
observed Gaia flux will appear over-luminous when fitting with a single white
dwarf model, resulting in a large radius and low white dwarf mass (McCleery et
al., 2020). Another scenario is when a companion strips the mass of a white
dwarf during a common-envelope evolution, leaving a white dwarf with a mass
that is smaller than what is possible assuming single star evolution (see, for
example, Temmink et al. 2020). A third path for binary evolution is a merger,
resulting in an apparently single white dwarf of any mass (Temmink et al.,
2020). We will address white dwarf mergers in a later section when describing
our simulation.

Unrelated to the above point, from Fig. 1 we can clearly see that there is a
downturn to lower masses for $T_{\rm eff}\lesssim 6000$ K or a cooling age of
about 2.5 Gyr. It is better illustrated by the purple solid line in Fig. 1,
which represents the median mass in a given $T_{\rm{eff}}$ bin. It is not
consistent with evolutionary models of white dwarfs, which predict cooling at
almost constant total mass and radius. In particular, Galactic population
simulations from Tremblay et al. (2016) have shown that the mean mass of white
dwarfs is predicted to stay the same within 1% for almost all look-back times,
except for the first $\approx$ 500 Myr after the Galactic disc formation,
where the mean mass of white dwarfs was slightly higher than it is now because
they came on average from slightly more massive progenitor stars.

The explanation for the small masses of cool white dwarfs is likely attributed
to an opacity which is missing or is incorrect in the atmospheric models
(Bergeron et al., 2019). In the absence of appropriate grids of model
atmospheres, we take the following ad-hoc approach. For each bin below 6000 K,
we find the difference between the median mass in the bin, and the overall
median mass above 6000 K, but below 20 000 K. We then apply this difference to
the masses of all white dwarfs below 6000 K as indicated by blood-orange
points in Fig. 1. This means that the median mass below 6000 K is then
corrected to agree with evolutionary models. This corrected median mass is
denoted in Fig. 1 by the dashed orange line. The specific equation for the
mass correction can be found in Tab. 1.

Now that white dwarf masses are corrected for the missing opacity, we can
remove any white dwarf whose corrected mass is still below the single star
evolution lower mass limit at 0.54 M⊙. This removes 110 white dwarfs, leaving
us with a sample of 973 white dwarfs.

Table 1: The mass correction applied to correct for issues with opacity. The table shows the amount of mass that needs to be added to individual white dwarfs for specific $T_{\rm{eff}}$ ranges, in order to get the median mass to agree to the median masses above 6000 K, but below 20 000 K, as expected from evolutionary models. Correction in mass |  $T_{\rm{eff}}$ range correction
---|---
(M⊙) | applies to (K)
$+$0.208 |  $3000\leq$ $T_{\rm{eff}}$ $<4000$
$+$0.138 |  $4000\leq$ $T_{\rm{eff}}$ $<5000$
$+$0.0395 |  $5000\leq$ $T_{\rm{eff}}$ $<6000$
Figure 1: The mass-$T_{\rm eff}$ distribution of the all-sky 40 pc white dwarf
sample. The black dots represent the original masses of the white dwarfs, with
the solid purple line representing the median mass in bins of around 1000 K.
Due to the unrealistic downturn in the median mass below 6000 K, the masses
had to be corrected. The corrected masses are denoted by blood orange dots
(see Table 1), with the corrected median mass indicated by the dashed orange
line.

In Fig. 2 we show a histogram of the absolute magnitude data we will be using
to determine the local stellar formation history. This is essentially the
white dwarf optical wavelength-range luminosity function for the 40 pc sample.
Given the precision of Gaia, the errors on absolute $G$ magnitude are
insignificant, with a median value of 0.03%. Instead, we plot the Poisson
errors. As expected, we see a sharp drop-off at dimmest magnitudes, which is
due to the finite age of the Galaxy. The histogram seems to have two peaks,
one centred around $G=12$ and the other around $G=14.5$. This has been seen in
other white dwarf luminosity functions (see, for example, Harris et al. 2006)
and is due to the different cooling mechanisms white dwarfs experience as they
age (Koester & Chanmugam, 1990).

Figure 2: A histogram of the absolute $G$ magnitudes for the 40 pc white dwarf
sample. The histogram is plotted in solid pink. Given the precision of Gaia
data, the individual errors on absolute $G$ magnitudes are insignificant.
Thus, we plot the Poisson errors in purple to account for the dominant source
of errors. The sizes of the bins are 0.5 dex in $G$ magnitude.

###  2.1 Separating the different Galactic populations

To derive the stellar formation history, we must first have some idea of the
shape it might take. As there are three Galactic components that can exist in
our local volume, the question arises whether they each have their own
separate history. To answer this, we must first see if the three components
are apparent in our sample. There are several methods used in literature to
identify and separate the populations of the thin disc, the thick disc, and
the halo of the Galaxy. They rely on the assumption that the components were
formed at different times, but ages can be difficult to measure for main-
sequence stars. Hence, the proxies of age, such as kinematics and/or
metallicity can be used to differentiate the populations, since, for example,
older stars are formed in lower metallicity environments and have larger
velocity dispersions than younger stars (Aumer & Binney, 2009; Casagrande et
al., 2011). When relying on observations of white dwarfs, the metallicity
information of the main sequence star is lost due to the high gravitational
fields of white dwarfs, which cause all but the lightest elements to sink out
of the atmosphere towards the core. Therefore, in the following we aim to use
the kinematics of white dwarfs to assign them to a specific Galactic
population. For the halo population, the typical Galactic velocities in the
local standard of rest are around 180-220 km s-1 (see, for example, Majewski
1992; Smith et al. 2009; Nissen & Schuster 2010; Du et al. 2018 and references
therein). The velocities of the thick disc stars are expected to be smaller
than the halo velocity, but larger than 70 km s-1 (Bensby et al., 2014). Thin
disc stars are considered to have velocities below 70 km s-1. Thus, we will be
using these values in the following analysis. In particular, we will use 180
km s-1 for the halo population in order to determine the maximum possible
number of halo white dwarfs.

From Gaia eDR3 we have access to the proper motions, right ascensions,
declinations, and distances for the white dwarfs in the 40 pc sample. If we
also have access to their radial velocities, we can derive their Galactic
velocity, where the component $U$ is directed towards the Galactic center,
component $V$ is in the direction of Galactic rotation, and component $W$ is
perpendicular to the Galactic disc. These components are corrected for the
local standard of rest. The first step is therefore to compile all available
radial velocities for our sample. To achieve this we used the radial velocity
catalogues of Anguiano et al. (2017) and Napiwotzki et al. (2020)333Note that
we also looked at samples of Pauli et al. (2006) and Raddi et al. (2022), but
in the former case, the radial velocities were older measurements and had
bigger errors than the ones found from Anguiano et al. (2017) and Napiwotzki
et al. (2020), and in the latter case there were no radial velocities for our
white dwarfs., which have been corrected for gravitational redshifts of white
dwarfs. In total, we found 64 (5.9% of total sample) northern and 59 (5.4% of
total sample) southern white dwarfs with radial velocities, all with hydrogen-
dominated atmospheres. Therefore, 11.4% of the white dwarfs in the 40 pc
sample had available radial velocities. If both Anguiano et al. (2017) and
Napiwotzki et al. (2020) samples had available radial velocities for a given
white dwarf, we chose the Napiwotzki et al. (2020) data, due to the sample’s
smaller errors and recentness. With the available radial velocities, we
calculate the Galactic velocities by utilising Astropy and its ICRS, SkyCoord,
GalacticLSR and Units packages and classes (Astropy Collaboration et al.,
2013, 2018). For white dwarfs with no available radial velocities, we instead
randomly sample a velocity from the radial velocity distribution of Raddi et
al. (2022). Once a radial velocity is assigned, we follow the same procedure
as before to find the Galactic velocity components. However, to reduce any
impact of randomness, we sample 1000 radial velocities for each white dwarf,
and use bootstrapping to determine the median values of the Galactic velocity
components and their associated errors. The errors are determined at the
confidence limits of 2.5% and 97.5%.

Fig. 3 shows the Toomre diagram for our 40 pc sample. The lower limits of the
total velocity for the thick disc and the halo populations are plotted as
dotted and dashed cyan lines, respectively. The Gaia errors are smaller than
the individual data points and the main source of error comes from the
observed radial velocity data, with the median error being around 8%. In terms
of the simulated radial velocities, the distribution of Raddi et al. (2022) is
well peaked, leading to a small error on the median Galactic velocities. To
test the effect of the more uncertain radial velocities, we also computed
Toomre diagrams where all the radial velocities were set to 0 km s-1 and where
the Raddi et al. (2022) distribution is shifted by several tens of km s-1. We
find no significant difference in the overall appearance of the plot, unless
the shift is of the order of 50 km s-1, which is improbable as it would make
the simulated velocities much higher than the velocities for the white dwarfs
with observed radial velocities.

Figure 3: The Toomre diagram for the 40 pc sample of white dwarfs. The plotted
velocity components are Galactic and have been calculated by taking into
account the local standard of rest. The dashed and dotted cyan lines denote
the total velocities equal to 180 km s-1 and 70 km s-1, which are reference
values that can be used to identify the halo and thick disc populations,
respectively. The colour of the individual points illustrates the derived
total age of a given white dwarf, as indicated by the right colour bar. The
green markers symbolize white dwarfs whose ages were larger than the age of
the Galaxy. Black contours represent white dwarfs with radial velocities found
from Anguiano et al. (2017) and Napiwotzki et al. (2020), while objects with
no contours have simulated radial velocity from the sample of Raddi et al.
(2022).

The colours of the markers in Fig. 3 indicate the total ages of the white
dwarfs. The white dwarf ages were derived using the University of Montréal
cooling tables (Bédard et al.,
2020)444http://www.astro.umontreal.ca/~bergeron/CoolingModels. In the 40 pc
sample, each white dwarf has a determined atmospheric composition from
observed spectra, therefore the appropriate cooling table (thick or thin
hydrogen layer) was used to determine the cooling age. Each white dwarf also
has a mass derived from Gaia data (corrected for opacity issues in the
previous section), which is transformed into the mass of its progenitor using
the initial-to-final mass relation of El-Badry et al. 2018. This progenitor
mass allows for the determination of the main sequence lifetime using the
Hurley et al. (2000) relation at solar metallicity, which is then combined
with the cooling age to find the total age. Some white dwarfs are denoted by
green colours as their ages could not be determined. This is because either
the atmospheric parameters were not available from the samples of McCleery et
al. (2020) and Tremblay et al. (2020); or their total age was larger than the
age of the universe. In total 78 white dwarfs do not have age determinations.

Overall, the Toomre diagram shows that all three potential regions of Galactic
components are populated, with the regions of the thin and the thick disc
components being most populous. This is expected because the thin disc is
concentrated in the Galactic plane where the Sun is found, while the thick
disc and the halo components are more spread out. The figure also indicates
that white dwarfs with lower velocities are more likely to be young. This is
expected because of the age-velocity dispersion relation. We must note that by
using the Toomre diagram we can identify potential candidates of the different
Galactic components, but we cannot separate them with absolute certainty. The
diagram seems to indicate there are two distinct velocity populations in our
sample: one concentrated ($\lesssim$100 km s-1) and one diffuse component
($\lesssim$300 km s-1) in the velocity space, with the diffuse component
mostly consisting of old stars (total age $\gtrsim$ 6 Gyr). From the data
available, it is difficult to assign these components to specific Galactic
populations such as thin and thick discs or halo, and these populations may
overlap in velocity space. As mentioned in the Introduction, other studies
have also found an overlap between the thin and thick discs in the local
neighbourhood, possibly due to the Gaia-Enceladus merger. Additionally, the
white dwarf study of Torres et al. (2019), which used Gaia data, has found 21
halo white dwarfs within 40 pc, suggesting the presence of the three Galactic
components in this volume. That is more halo white dwarfs than we find,
however, their classification technique is more advanced than the technique
used in this paper. In some cases they classify white dwarfs with velocities
smaller than 180 km s-1 as halo white dwarfs, based on their other intrinsic
parameters. We therefore proceed by testing stellar formation histories with
one, two and three Galactic components to account for all possible scenarios.

###  2.2 Age-velocity dispersion

The age-velocity dispersion relation introduces an observational bias which we
must consider when comparing Galactic simulations with the observations. Older
white dwarfs are statistically more likely to have larger velocities and as
such are less likely to be observed in a given volume close to the Galactic
plane. If this bias was not taken into account, it could lead to an erroneous
conclusion that the stellar formation history must have peaked at recent times
to explain an apparent over-abundance of younger white dwarfs. In this work we
consider this bias to affect the $W$ component of the total velocity only, as
we assume to first order that the same amount of objects come in to and out of
40 pc volume in the directions of the $U$ and $V$ velocity components due to
radial mixing.

In Fig. 4 we plot binned values of velocity dispersion and total age for the
40 pc white dwarf sample. We use ten bins with nearly equal number of white
dwarfs. We remove white dwarfs whose ages could not be determined. For each
bin we find the median total age and its error using bootstrapping, where the
value of the error is taken as the confidence intervals of 95.7% and 2.5%. To
determine the velocity dispersion, we also use bootstrapping with the same
confidence intervals as errors. Bootstrapping is beneficial in our case, as
the errors on individual measurements are almost insignificantly small when
compared to the actual scatter of the data points. Bootstrapping also allows
us to easily account for the random fluctuations arising from using simulated
radial velocities.

Figure 4: The binned age-velocity dispersions for the 40 pc sample of white
dwarfs. These relations are plotted for the three Galactic velocity
components, $U$, $V$ and $W$, in solid orange, dashed pink, and dotted black
lines, respectively. Median total ages and errors are found for each bin using
bootstrapping. The velocity dispersions are also found using bootstrapping,
taking into account that for most white dwarfs the radial velocities were
simulated using the Raddi et al. (2022) radial velocity distribution. The
errors shown on the plot represent the confidence intervals of 2.5% and 97.5%.

In line with other studies (Tremblay et al., 2014; Tremblay et al., 2016;
McCleery et al., 2020; Raddi et al., 2022), Fig. 4 shows the velocity
dispersion increasing with age. The range of the velocity dispersion, 10-60 km
s-1 is also consistent with the Raddi et al. (2022) sample. Additionally, we
find that the $U$ component has the largest velocity dispersion. Overall, we
find that the velocities of the 40 pc white dwarf sample resemble previous
studies of white dwarfs, and the main sequence stars, but further study is
needed, since not all white dwarfs in our sample have observed radial
velocities. We include the age-velocity dispersion in our simulation, as
described in the next section, Sect. 3.

##  3 Numerical simulation of the 40 pc white dwarf sample

The simulation of the local white dwarf sample is based on Monte Carlo
methods. To ensure that the random number fluctuations are insignificant, we
must generate a large enough number of simulated white dwarfs, which we fix at
3000 white dwarfs. We find that the fluctuations between the different
simulations generated with this number of white dwarfs are on 2% level. We
define the total age of the star as the time difference from the present time
and the time of stellar formation (formation time). The time variable used
throughout this work is the lookback time, going back from present day ($t=0$)
to formation time. In the simulation, the formation time is assigned based on
the assumed stellar formation history. To generalise the stellar formation
history in a way that is suitable for fitting observed data, we assume either
one, two or three Galactic components that each have their own uniform stellar
formation history. The specific details on this can be found in Sect. 4.

For each white dwarf, a random position is generated in the three Cartesian
($x,y,z$) coordinates. Both the $x$\- and $y$-coordinates are sampled from a
uniform distribution which spans the distances between 0 and 40 pc from the
Sun. As mentioned in Sect. 2.2, we assume to first order that the amount of
stars coming in and out of the 40 pc in $x$ and $y$ directions are equal. On
the other hand, the $z$-coordinate is sampled from a quasi-exponential
distribution, which has been offset by 20 pc, due to the Sun’s vertical
position relative to the Galactic centre (see, for example, Siegert 2019),
i.e. from the probability density

$\rho{(z)}=1-\exp{(-\mathrm{(z+20})/\mathrm{scale\ height)}}.$ (1)

The $z$-distribution is limited to distances between 0 and 40 pc. An important
assumption in this distribution is the scale height of the Galactic component.
In Sect. 2.2 we found the relation between the velocity dispersion of the $W$
velocity component as a function of total age of the star. Assuming that at
the age of 1 Gyr the scale height is 75 pc (Wegg & Phinney, 2012; Buckner &
Froebrich, 2014) and a simple monotonic function fit to Fig. 4 (green curve),
we can relate the total age of the star directly to the scale height using the
$W$ velocity dispersion relation. We show our derived scale height as a
function of time in Fig. 5. Thus, the older the star is, the higher its scale
height, and as such, older white dwarfs are less likely to be found within the
40 pc volume for the same stellar formation history. Note that, our velocity
dispersion relation for the $W$ velocity component only goes up to around 7
Gyr. If the formation time is larger than this value, we assume that the scale
height is the same as the scale height at the last available age of the
observed velocity dispersion relation. Once the $x,y$\- and $z$-coordinates
are generated, any star that is more than 40 pc from the Sun is rejected, and
the process is repeated until a star within 40 pc is created.

Figure 5: The scale height as a function of total age of the star. This
relation has been derived based on a smooth monotonic fit of observed velocity
dispersion of the $W$-velocity component, and the assumption that at age 1 Gyr
the stars have vertical scale height of 75 pc (Wegg & Phinney, 2012; Buckner &
Froebrich, 2014).

An initial mass is then assigned for each star, assuming a mass distribution
proportional to $M_{*}^{-2.3}$, where $M_{*}$ is the initial mass (Kroupa,
2001, 2002). By using the main sequence age relation of Hurley et al. (2000),
the main sequence lifetime can be derived from its initial mass and from the
assumed solar metallicity, which we set as $Z=0.0134$ (Asplund et al., 2009).
Note that this is the present-day photospheric metallicity of the Sun and not
the "proto-Solar" interior metallicity at the moment of the formation of the
star. Through the initial-to-final mass relation, IFMR, the initial mass of
the star gives us the final white dwarf mass. In the default simulation we use
the IMFR of El-Badry et al. (2018). We will address these assumptions further
in Sect. 6.

The cooling age of the white dwarf is defined as the difference between the
formation time of the initial star (total age) and the main sequence lifetime.
If the resultant white dwarf cooling age is negative (including any merger
delay or hastening as discussed in Sect. 3.1), it means that there has not
been enough time for the star to evolve past the main-sequence and form a
white dwarf. Therefore, the process is restarted until a suitable white dwarf
candidate is found.

The mass, cooling age and absolute $G$ magnitude of the simulated white dwarf
is based on evolutionary models from Bédard et al. (2020) and appropriate
model atmospheres with H- or He-composition (Tremblay et al., 2011;
Cukanovaite et al., 2021). In our simulation we assume that 25% of all white
dwarfs have helium-dominated atmospheres, hence thin hydrogen envelopes and
appropriate synthetic magnitudes for DB or DC white dwarfs, which is based on
the observed number from the 40 pc sample (McCleery et al., 2020). Other
studies have found similar numbers (see, for example, López-Sanjuan et al.
2022). At this point, the simulation has produced an absolute $G$ magnitude
distribution that can be compared against real data.

###  3.1 Binary system interactions

In the simulation, we also take into account time delays which arise from
binary mergers. To our knowledge, previous studies of stellar formation
history involving white dwarfs did not address this delay. To take this into
account in our simulation, we use the results of the binary population
simulations of Temmink et al. (2020), namely the merger progenitor fractions
and the time delays or speed ups due to the mergers. In particular we use the
default models from Temmink et al. (2020) Figs. 6 and 9, which we replot in
Fig. 6. Temmink et al. (2020) employ the binary population synthesis (BPS)
code SeBa (Portegies Zwart & Verbunt, 1996; Toonen et al., 2012) to simulate
large numbers of single stars and binary systems. In SeBa, processes such as
stellar winds, mass transfer, angular momentum loss, common envelope phases,
gravitational radiation and stellar mergers are considered with the
appropriate prescriptions. Details of the setup of the default model can be
found in Temmink et al. (2020), but we briefly summarize the most relevant
parts in the following. The initial binary fraction is assumed to be 50%, and
the authors model common envelope phases with the alpha-prescription (with
$\alpha\times\lambda=2$), which is based on the energy budget in the binaries.
The authors then consider all single stars and binary systems that eventually
form a single white dwarf. The most important contributing channels can be
seen in Figs. 1 and 3 of Temmink et al. (2020). Here, we use the distributions
they find in populations formed after a burst in star formation.

Figure 6: Top figure: The fraction of white dwarfs that have formed via a
merger as a function of the mass of the white dwarf. This is the default model
from Temmink et al. (2020). Bottom figure: The time difference in Gyr between
the prior history of a white dwarf assuming it has formed via binary
evolution, $\tau_{\mathrm{BE}}$, and single star evolution,
$\tau_{\mathrm{SSE}}$, as a function of white dwarf mass. The prior history is
the time taken to form the final observed white dwarf, all the way from the
formation of the original star, including any subsequent merger events in the
case of binary evolution model. The pink horizontal lines indicate the median
value of the time difference. The black boxes indicate the first and third
quartiles of the time difference distributions. The values shown here are from
the default model of Temmink et al. (2020).

If any of the apparently single 40 pc white dwarfs are formed via a binary
merger, their prior history is different to the history of a white dwarf
formed through single star evolution. As shown in the bottom plot of Fig. 6,
in most cases this will mean that its prior history takes a longer time, but
for lower mass white dwarfs it can actually be faster. Before we can apply
this time correction, we first need to know how many white dwarfs in our
simulation have formed as a result of a merger. This fraction is shown in the
top panel of Fig. 6. It is apparent that more massive white dwarfs are more
likely to be formed as a result of a merger, but even at lower masses such as
0.6 M⊙, the fraction is as high as 20% for the default model of Temmink et al.
(2020). Therefore, Fig. 6 (top panel) gives the probability of a simulated
white dwarf with a given mass being the product of a merger, in which case we
apply a time delay based on the delay distribution shown in the bottom panel
of Fig. 6. We determine the time correction by randomly sampling from a
Gaussian distribution with a mean and standard deviation determined from the
quartile information (Q1, Q2 and Q3) given in Fig. 6 of Temmink et al. (2020).
Note that this corresponds to the combination of a cooling delay and extended
pre-white dwarf lifetime. Since the time difference is almost always positive,
in most cases it will shorten the cooling age of the white dwarf for the same
stellar formation history. Therefore, when compared to a model with no binary
evolution, the binary-corrected model should have hotter and brighter white
dwarfs for the same formation history. Binary evolution is included in our
simulation by default and is used in the following analysis.

We remind the reader that we also account for binarity by removing the
majority of double degenerate candidates from our mass cut at 0.54 M⊙
introduced in Section 2. Unresolved white dwarf and main-sequence pairs are
also missing from our sample as it was constructed from a Gaia selection of
single or multiple white dwarfs (Gentile Fusillo et al., 2021).

##  4 Determination of Stellar Formation History

We fit our simulated model, with variable stellar formation history, to the
observed absolute $G$ magnitude distribution of the 40 pc white dwarf sample.
We choose our fitting data to be the empirical cumulative distribution
function, ECDF, of absolute $G$ magnitude, which unlike the histogram shown in
Fig. 2 does not depend on arbitrarily chosen bin number or size. Instead, the
number of bins is determined by the number of observed white dwarfs. The ECDF
of the 40 pc white dwarf sample is shown in Fig. 7. The errors on the absolute
$G$ magnitude are small enough to be invisible in the figure, therefore, we
use propagated Poisson errors.

Figure 7: Top plot: The empirical cumulative distribution function of the
absolute magnitudes. The lime green error bars show the 40 pc sample ECDF,
whereas in solid black, dashed purple, and dotted orange we plot the best fit
ECDFs found by assuming uniform stellar formation histories of either one, two
or three Galactic components. Middle plot: The residuals from the fits
compared to the observational errors. The errors are indicated as filled green
area. Bottom left plot: The best fitting stellar formation history assuming
one Galactic component. In solid black we show all the stellar formation
histories that correspond to the probability distribution function of the
fitting parameter, $\tau_{\mathrm{max}}$. In dashed orange we plot the stellar
formation history found from the median value of the probability distribution
function of $\tau_{\mathrm{max}}$, which we define as the best fit stellar
formation history found from the observed 40 pc sample. Bottom centre plot:
This plot is similar to the bottom left plot, but shows the best fit two
Galactic component stellar formation history. Bottom right plot: Same as
bottom left and centre plots, but this plots shows the best fit three Galactic
component stellar formation history.

For a given simulation, i.e. given stellar formation history, we calculate the
ECDF using all 3000 simulated white dwarfs. To directly compare our simulation
to the data, we have to scale the simulated ECDF to the number of white dwarfs
in our data sample, i.e., $\sim$1000\. We do this by interpolating the
simulated ECDF at the values of observed 40 pc absolute $G$ magnitudes. As
mentioned previously, each simulation depends on a stellar formation history,
which we parameterise in two ways. In the first case, we assume that there is
only one Galactic component in the simulation, and that the stellar formation
history is uniform between present day (formation time of 0 Gyr) and some
maximum time, $\tau_{\mathrm{max}}$. When fitting this scenario to the
observed 40 pc sample, we only fit for the parameter $\tau_{\mathrm{max}}$.

In the two Galactic component case, we assume that an older component starts
to uniformly produce stars at a time, $\tau_{\mathrm{max,2}}$ until the time,
$\tau_{\mathrm{max,1}}$. At this point the older component shuts down, and the
younger components starts producing stars from the time
$\tau_{\mathrm{max,1}}$ until present day. In this case, we also have a
parameter $N_{2}$, which describes the fraction of stars in the older
component. Note that the fraction of stars in younger component will be
1-$N_{2}$. Therefore, in this case we have three fitting parameters:
$N_{2},\tau_{\mathrm{max,1}},\tau_{\mathrm{max,2}}$. The three Galactic
component case is similar, but the fitting parameters are
$N^{{}^{\prime}}_{2},N^{{}^{\prime}}_{3},\tau^{{}^{\prime}}_{\mathrm{max,1}},\tau^{{}^{\prime}}_{\mathrm{max,2}},\tau^{{}^{\prime}}_{\mathrm{max,3}}$,
where the subscript number refers to either the first, second or third
Galactic component.

As an aside, note that the absolute stellar formation rate is set by the total
number of white dwarfs within 40 pc, and so is not a free parameter. Since the
vast majority of local stars are M dwarfs that have not yet formed white
dwarfs, we refrain from the derivation of the absolute stellar formation rate
per unit volume and time. Furthermore, we do not correct for the larger
velocity of halo stars, and as such we make no attempt at deriving the
absolute stellar halo formation rate. This would be difficult to constrain in
any case since the halo white dwarf fraction within 40 pc is smaller than 2.5%
(Torres et al., 2019).

In each case, we fit the data by minimising the reduced weighted $\chi^{2}$,

$~{}\chi^{2}=\sum_{G^{\prime}=\min(G_{\mathrm{observed}})}^{G^{\prime}=\max(G_{\mathrm{observed}})}\frac{\mathrm{(simulated\
ECDF}_{G^{\prime}}-\mathrm{observed\ ECDF}_{G^{\prime}})^{2}}{\mathrm{error\
on\ observed\ ECDF}_{G^{\prime}}^{2}\times\mathrm{dof}},$ (2)

where $G^{\prime}$ is the index for a given absolute magnitude bin of the
ECDF; and dof is degrees of freedom, defined as

$\mathrm{dof}=n_{\mathrm{data}}-(n_{\mathrm{parameters}}-1),$ (3)

where $n_{X}$ is either the number of data points or number of fitted
parameters.

For any fitting procedure, the first step is to determine the initial
parameters. To do this, we explore the entirety of the physical parameter
space by calculating $\chi^{2}$ for select number of stellar formation
histories. The parameters that give the smallest $\chi^{2}$ are then chosen as
initial parameters for the minimisation. The minimisation is based on the
Scipy minimize package in Python (Virtanen et al., 2020; Van Rossum & Drake,
2009) and is combined with our custom MCMC program to find the probability
distributions of the best-fitting parameters. The final best fit parameters
are then taken as the median of the probability distribution with the errors
quoted being the 2.5% and 97.5% confidence limits. When assuming a one-
component stellar formation history the best fit parameter is found to be
$\tau_{\mathrm{max}}=10.6^{+0.5}_{-0.5}$ Gyr. For the two component stellar
formation history the best fit parameters are $N_{2}=0.39^{+0.03}_{-0.04}$,
$\tau_{\mathrm{max,1}}=6.6^{+0.5}_{-0.4}$ Gyr, and
$\tau_{\mathrm{max,2}}=10.4^{+0.7}_{-1.1}$ Gyr. For the three Galactic
component case, the best fit parameters are
$N^{{}^{\prime}}_{2}=0.29^{+0.02}_{-0.02},N^{{}^{\prime}}_{3}=0.25^{+0.01}_{-0.02},\tau^{{}^{\prime}}_{\mathrm{max,1}}=5.0^{+0.3}_{-0.3},\tau^{{}^{\prime}}_{\mathrm{max,2}}=8.4^{+0.5}_{-0.6},\tau^{{}^{\prime}}_{\mathrm{max,3}}=11.2^{+0.7}_{-0.7}$.

The best fits are shown visually in Fig. 7. The top panel shows the observed
and model ECDFs, the middle plot shows the residuals of the fits compared to
the observational errors, and the bottom plot shows the corresponding SFHs for
the one, two or three Galactic component case. Note that since we use MCMC to
derive the best fit parameters, we find their probability density functions.
These can be converted to the probability density functions of the best fit
stellar formation histories. These are shown as black lines in the bottom plot
of Fig. 7, where the orange dashed lines indicate the median value of the SFH,
which we take to be our best fit SFH. While looking at the top and middle
plots of Fig. 7 we can see that there are no significant differences between
the fits. Even the best fit SFHs shown in the bottom are all similar, with
only minor deviations from complete uniformity for the SFHs with two or three
Galactic components. Beyond the work shown in this paper, we have tried fits
with skewed Gaussians instead of uniform distributions, both with one, two and
three Galactic components. These more complex stellar formation histories do
provide fits on par or in slightly better agreements with the 40 pc
observations. This is not surprising, since more complex models with larger
number of parameters allow for more molding of resultant ECDF. As such, better
fitting to the data could be due to over-fitting. Therefore, our philosophy is
to assume the simplest solution which gives a reasonable fit to the data, i.e.
a uniform SFH with one Galactic component. We will show this in the following
section by comparing it to other stellar formation histories from literature.

##  5 Comparison to other stellar formation histories

In the left plot of Fig. 8 we compare the ECDFs found from our simulations
with different stellar formation histories. These histories were taken from
Reid et al. (2007), Tremblay et al. (2014), Fantin et al. (2019) and Mor et
al. (2019) studies, and are shown in the right plot of Fig. 8. The stellar
formation histories of Reid et al. (2007) and Mor et al. (2019) are based on
main sequence stars, whereas Tremblay et al. (2014) and Fantin et al. (2019)
are based on white dwarfs. The Mor et al. (2019) and Fantin et al. (2019)
results are based on Gaia data. We chose these particular SFHs as they have
different shapes, specifically they show peaks at different times.

Figure 8:  Top plot: A comparison between our one Galactic component stellar
formation history and stellar formation histories found from literature. The
SFH of Reid et al. (2007), Tremblay et al. (2014), Fantin et al. (2019) and
Mor et al. (2019) are plotted in dashed-dotted purple, dashed blue, dotted
pink and solid orange, respectively. Our SFH is plotted in solid black. All
stellar formation histories are normalised. Middle plot: The observed ECDF
plotted as green errorbars is compared to simulated ECDFs when the SFH in our
simulation is changed to the SFHs of Reid et al. (2007), Tremblay et al.
(2014), Fantin et al. (2019) and Mor et al. (2019). These are plotted in solid
black, dash-dotted purple, dashed blue, dotted pink and solid orange,
respectively. Bottom plot: The residuals of the above plot.

From this comparison we find that our uniform SFH and the SFHs from Reid et
al. (2007) and Mor et al. (2019) fit the data, whereas the SFHs from Tremblay
et al. (2014) and Fantin et al. (2019) do not agree with the data and can be
excluded. The data from the 40 pc sample favours a more flat SFH such as ours
or Reid et al. (2007), or a SFH with two peaks such as the SFH of Mor et al.
(2019).

In Fig. 9 we show a comparison of the absolute $G$ magnitude histograms from
the observed 40 pc white dwarf sample and simulations with published stellar
formation histories of Reid et al. (2007), Tremblay et al. (2014), Fantin et
al. (2019) and Mor et al. (2019). The residuals of the comparisons are shown
in the bottom plots of Fig. 9. As for the case of the empirical cumulative
distribution function of absolute $G$ magnitudes, the best agreement is with
our own stellar formation history, as well as those of Reid et al. (2007) and
Mor et al. (2019).

Figure 9: Top left plot: A histogram of the absolute $G$ magnitudes for the 40
pc sample is plotted in solid pink. The histogram is normalised to unit area.
The absolute $G$ magnitude histogram found from the best fit of our one-
component simulation is plotted in yellow. In green and blue we plot the
absolute $G$ magnitudes from our simulation relying on the stellar formation
history of Tremblay et al. (2014) and Fantin et al. (2019), respectively. The
residuals of this comparison are shown in the bottom left plot, with the solid
pink line indicating perfect agreement. Top right plot: Same as the top left
plot, but in yellow and blue we plot the histograms found from our simulation
with the stellar formation histories taken from Reid et al. (2007) and Mor et
al. (2019). The bottom right plot shows the corresponding residuals.

In the following, we further test the robustness of our stellar formation
history by assessing how various assumptions made in our models affect the
best fit.

##  6 Theoretical assumptions

Assumptions have been made when building our simulation. These include complex
assumptions based on difficult physical considerations, such as the relation
between the initial and final mass of a star as it becomes a white dwarf. In
the following, our aim is to assess whether the assumptions made in our
simulations have a significant effect on the derived stellar formation
history. We do this by changing one assumption at a time and assessing how
this change affects the best-fit $\chi^{2}$ between data and simulations. Note
that the value of $\chi^{2}$ can indicate whether a given hypothesis is in
agreement with the data. For this analysis we proceed with a one-component,
uniform stellar formation history as we have shown that it gives a reasonable
fit to the data.

###  6.1 Metallicity

In the simulation, the assumed metallicity has a direct effect on the derived
main sequence lifetime, which depends both on the metallicity and on the main
sequence mass. In general, for main sequence masses below $\approx 6$ M⊙,
metallicities smaller than the solar metallicity result in shorter main
sequence lifetimes. If the formation time does not change, i.e. we use the
same stellar formation history, then the resulting white dwarfs will have
longer cooling ages and dimmer magnitudes. At around 6 M⊙ and above, the
metallicity has the opposite effect on the main sequence lifetime. Due to the
initial mass function, most main sequence stars in our simulation will have
lower masses, i.e., below 6 M⊙. Thus, the reversal of the effect of
metallicity above 6 M⊙ is not significant for our simulation.

When trying to fit data with our simulation, the effect of metallicity on our
SFH and therefore the derived value of the parameter $\tau_{\mathrm{max}}$ is
the following: A simulation with larger assumed metallicity will produce more
bright white dwarfs than a simulation with low metallicity if the stellar
formation history stays the same. Therefore, to compensate for this, the
fitted value of $\tau_{\mathrm{max}}$ should be larger for the simulation with
larger metallicity in order to produce enough dim white dwarfs. Similarly, it
could change the shape of the fitted SFH from uniform to a peak at later
times, to over-compensate for too many simulated bright white dwarfs. We run
three additional types of simulations with different assumed metallicities:
$0.2+Z_{\odot}$, $Z_{\odot}/10$ and $Z_{\odot}/1000$, where $Z_{\odot}$ is the
solar metallicity equal to 0.0134, used in the original simulation. These
particular choices of metallicity are extreme and have been chosen to test the
robustness of our derived stellar formation history. The metallicity of
$0.2+Z_{\odot}$ represents metal-rich stars, which have metallicity around 0.2
dex higher than the solar metallicity (Feltzing & Gonzalez, 2001). The
metallicities of $Z_{\odot}/10$ and $Z_{\odot}/1000$ represent the range of
the metallicities found in metal-poor population II stars. Note that for this
paper we assume that the metallicity is constant within the 40 pc sample, as
Rebassa-Mansergas et al. (2021) showed from white dwarfs with wide stellar
companions and within 500 pc of the Sun that there is no evidence of
metallicity changing as a function of formation time.

The 40 pc sample ECDF is fitted with the three types of simulations using the
MCMC method with 1000 walkers. Each walker finds a best fit value of
$\tau_{\mathrm{max}}$ and the corresponding $\chi^{2}$ value. Thus, for each
simulation we find a probability density function of both $\chi^{2}$ and
$\tau_{\mathrm{max}}$. In Table 2 we list the median values of $\chi^{2}$ and
$\tau_{\mathrm{max}}$, and as their errors we show the confidence limits of
2.5% and 97.5% of their distributions. As expected, the increased metallicity
results in a larger value of $\tau_{\mathrm{max}}$ when compared to other
metallicities and the results from the original simulation. The $\chi^{2}$
values stay close to 1 in all cases and thus indicate that different values of
metallicity still allow the uniform SFH to agree with the data, without
needing a more peaked SFH.

Table 2: The best fit $\chi^{2}$ and $\tau_{\mathrm{max}}$ found from simulations where a specific assumption is changed. The first column indicates which assumption is changed, with Metallicity, IMF, $\tau_{\mathrm{MS}}$, IFMR, Velocity correction, Binary correction, He-fraction, WD cooling model referring to tests with changes in either the metallicity, initial mass function, main sequence lifetime, initial-to-final mass relation, the velocity correction, binary correction, the fraction of helium-dominated atmosphere white dwarfs, the white cooling model, or the 40 pc sample from mass correction and double degenerate removal, respectively. The value of the assumption is listed in the second column, with more information provided in the text. Third column lists the best fit $\chi^{2}$ found from 1000 walkers, with the errors indicating the confidence limits of 2.5% and 97.5% of $\chi^{2}$ probability density functions. The last column is the same but for the parameter $\tau_{\mathrm{max}}$, which uniquely defines the uniform, one Galactic component SFH. Simulation | Assumption | $\chi^{2}$ |  $\tau_{\mathrm{max}}$ (Gyr)
---|---|---|---
Original |  | $1.4^{+0.8}_{-0.7}$ | $10.6^{+0.5}_{-0.5}$
Metallicity | $Z_{\odot}+0.2$ | $1.1^{+0.8}_{-0.5}$ | $11.0^{+0.5}_{-0.7}$
| $Z_{\odot}/10$ | $1.3^{+0.8}_{-0.6}$ | $10.0^{+0.4}_{-0.4}$
| $Z_{\odot}/1000$ | $1.5^{+0.9}_{-0.6}$ | $9.6^{+0.4}_{-0.4}$
IMF | Eq. 4 (steeper) | $2.6^{+2.0}_{-1.0}$ | $11.7^{+0.7}_{-1.1}$
| Eq. 5 (shallower) | $0.6^{+0.5}_{-0.3}$ | $9.6^{+0.4}_{-0.4}$
$\tau_{\mathrm{MS}}$ | $(1-0.048)*\tau_{\mathrm{MS}}$ | $1.4^{+0.9}_{-0.6}$ | $10.5^{+0.6}_{-0.2}$
| $(1+0.048)*\tau_{\mathrm{MS}}$ | $1.4^{+0.9}_{-0.6}$ | $10.6^{+0.6}_{-0.5}$
IFMR | Cummings et al. (2018) | $1.0^{+0.9}_{-0.5}$ | $11.1^{+0.6}_{-0.6}$
| Andrews et al. (2015) | $1.5^{+0.9}_{-0.7}$ | $9.6^{+0.5}_{-0.4}$
Velocity correction | Not included | $1.5^{+0.9}_{-0.7}$ | $10.6^{+0.6}_{-0.5}$
Binary correction | Not included | $1.5^{+0.9}_{-0.7}$ | $10.4^{+0.5}_{-0.6}$
He-fraction | 0% | $1.7^{+0.9}_{-0.7}$ | $10.6^{+0.5}_{-0.5}$
| 50% | $1.4^{+0.9}_{-0.7}$ | $10.6^{+0.5}_{-0.5}$
WD cooling model | LPCODE | $4.6^{+1.6}_{-1.2}$ | $9.7^{+0.6}_{-0.6}$
40 pc mass correction | No correction | $1.6^{+0.8}_{-0.6}$ | $8.6^{+0.5}_{-0.5}$
and double degenerate removal | No correction and no removal | $1.1^{+0.8}_{-0.5}$ | $10.3^{+0.5}_{-0.5}$

In the top left plot of Fig. 10 we show the residuals of the fits and how they
compare to the observed propagated Poisson errors, which are plotted as filled
green area. It is clear that all simulations give reasonable fits, but the
median value of $\chi^{2}$ is smallest for the simulation with
$Z_{\odot}+0.02$. However, looking at the top left plot of Fig. 10 it is clear
that the higher metallicity does not describe dimmer magnitudes as well as
other simulations, especially the $Z_{\odot}/1000$ simulation.

Figure 10: Left plots: The residuals from the fits of different simulations to
the 40 pc data. In all left plots the solid black line indicates the residuals
between the original simulation and the data. The green filled area denotes
the propagated Poisson errors. Top, centre and bottom left plots show the
residuals when the metallicity, the IMF, and the main sequence lifetimes are
changed, respectively. The legends indicate which coloured lines represent a
given assumption and its chosen value. Right plots: The one-component, uniform
stellar formation histories found from fitting the 40 pc sample with
simulations where a single assumption is changed. Each right plot corresponds
to the left plot on the same row. The right plots shows the full stellar
formation history, where the formation time of zero Gyr means present day. The
legends indicate the best fit stellar formation histories and which assumption
they correspond to. The corresponding coloured areas represent the 2.5% and
97.5% confidence limits of the derived stellar formation history when the
relevant assumption is changed. For clearer illustration of the differences in
the derived stellar formation history, the inset plot shows the zoomed in area
near the best fit formation time of $\tau_{\mathrm{max}}$.

In the top right plot of Fig. 10 we show the best fit stellar formation
history for each type of simulation, alongside the confidence limits of 2.5%
and 97.5%, which are shown as correspondingly coloured areas. There is a
significant overlap between the original simulation and the other three types
of simulations, suggesting that the metallicity has a small effect on the
derived stellar formation history within the framework of one Galactic
component uniform SFH. If we take into account that the chosen metallicities
are extreme choices, this conclusion becomes more robust, as it is unlikely
that within 40 pc the metallicity could change so drastically from solar
metallicity. We conclude that the assumed metallicity does not change our
conclusion that a one-component Galactic uniform SFH fits the 40 pc data well.

###  6.2 Initial mass function

In the original simulation we use the initial mass function (IMF) from Kroupa
(2001), namely, $\rho(M)\propto M^{-2.3}$, where $M$ is the initial mass of
the star and $\rho(M)$ is the initial mass density function. Kroupa (2001)
provides errors for the power index of the function, which depends on the mass
of the star. In Eqs. 4 and 5 we show the resultant initial mass density
functions when the errors are either added to or subtracted from the power
index of the Kroupa (2001) IMF.

$~{}\begin{gathered}\rho(M)\propto M^{-2.6}\mathrm{\ for\ }M<1M_{\odot},\\\
\rho(M)\propto M^{-3.0}\mathrm{\ for\ }M\geq 1M_{\odot},\end{gathered}$ (4)

$~{}\begin{gathered}\rho(M)\propto M^{-2.0}\mathrm{\ for\ }M<1M_{\odot},\\\
\rho(M)\propto M^{-1.6}\mathrm{\ for\ }M\geq 1M_{\odot}.\end{gathered}$ (5)

Eq. 4 is steeper than the original IMF used, whereas Eq. 5 is shallower. To
test whether the IMF affects the accuracy of the derived stellar formation
history, we run two additional simulations, where we replace the Kroupa (2001)
IMF with the Eqs. 4 and 5. Steeper IMF means that more low mass main sequence
stars will be created. Therefore, main sequence lifetimes will be on average
longer than in the original simulation, and if we assume the same stellar
formation history and thus formation times, the resultant white dwarf cooling
ages will be shorter. This means that for the same stellar formation history,
the simulated ECDF will have more white dwarfs at brighter magnitudes. This is
a similar effect to increasing metallicity as described in Sect. 6.1. The
opposite is true for the shallower IMF from Eq. 5.

In Table 2 we show the best fit values of $\chi^{2}$ and
$\tau_{\mathrm{max}}$. The much larger value of $\chi^{2}$ for the simulation
with steeper IMF indicates that for this assumption the uniform SFH does not
agree with the data within one standard error. For the shallower IMF, we get a
better fit than compared to the original model. These two conclusions are
supported by the central row left plot of Fig. 10. The central right plot of
Fig. 10 agrees with our expectations of steeper IMF producing more bright
white dwarfs and therefore resulting in larger $\tau_{\mathrm{max}}$ to
produce enough dim white dwarfs to match observations when compared to the
original simulation. Both figures show that the chosen IMF has a moderate
effect on the derived stellar formation history.

###  6.3 Main sequence age function

In the original simulation, the main sequence lifetime is found from the age
function of Hurley et al. (2000) which has an error of 4.8%. To test whether
the main sequence age function has a significant effect on the derived stellar
formation history, we run two additional simulations where the error of 4.8%
is either added to or subtracted from the age given by the Hurley et al.
(2000) relation. When compared to the original simulation, an increase in the
main sequence lifetime will result in shorter white dwarf cooling ages, since
the formation time remains the same. Therefore, such a simulation will have
more white dwarfs at brighter magnitudes in the resultant ECDF. Again, this is
similar to increasing the metallicity (Sect. 6.1) or having a steeper IMF
(Sect. 6.2). The opposite is expected for the simulation where the error is
taken away from the main sequence lifetime.

Table 2 shows the best fit $\chi^{2}$ and $\tau_{\mathrm{max}}$ when changing
the main sequence lifetime. Additionally, the bottom row of Fig. 10 visually
illustrates the results of this test. The best fit values of $\chi^{2}$ and
$\tau_{\mathrm{max}}$ with either test agree with the original simulation. The
bottom left plot shows that visually the simulated ECDFs from the three types
of simulations look identical. The probability distributions shown in centre
left plot are also almost identical. It is clear that the effect on the
resultant ECDF is insignificant, therefore, we proceed with using the original
main sequence lifetime relation from Hurley et al. (2000), and do not
investigate this assumption any further.

###  6.4 Initial to final mass relation

To test whether the IFMR has a significant effect on the derived stellar
formation history we run two additional simulations which replace the El-Badry
et al. (2018) IFMR with the Andrews et al. (2015) or Cummings et al. (2018)
relations, respectively. We choose these two particular IFMRs because they are
all relatively recent, but also give significantly different final masses of
white dwarfs. In most cases, the Cummings et al. (2018) relation results in
higher white dwarf masses than the El-Badry et al. (2018) relation, whereas
the Andrews et al. (2015) relations always gives smaller white dwarf masses
than the El-Badry et al. (2018) relation. Note that the simulated white dwarf
cooling age is not directly changed by the IFMR, since the main sequence
lifetime and the stellar formation history remain the same. As per white dwarf
cooling tables and model atmospheres, the white dwarf mass has a direct effect
on the Gaia $G$ magnitude. In general, an increase in white dwarf mass results
in dimmer magnitudes for all white dwarf cooling ages below $13$ Gyr.
Therefore, the simulation with Cummings et al. (2018) IFMR should produce more
dim white dwarfs than the simulation with El-Badry et al. (2018) IFMR. This is
similar to what we observed with previous assumptions.

Table 2 and the top row of Fig. 11 show the results of these tests. The
$\chi^{2}$ between the different types of simulations does not change
significantly, therefore all three simulations are broadly consistent with the
observed ECDF. The simulation with Cummings et al. (2018) IFMR results in
larger $\tau_{\mathrm{max}}$ than the original and the Andrews et al. (2015)
simulations. We note that the effect of the IFMR is similar to what we
observed with changing the initial mass function, except it is less
significant. Both these assumptions affect the simulated masses of white
dwarfs, hence these two parameters may be degenerate in determining the
stellar formation history.

Figure 11: This figure is similar to Fig. 10, but here we plot the results of
the tests where the initial to final mass relation is changed (top row); where
the velocity correction is removed (centre row); where the binary correction
is removed (centre row); and where the fraction of He-atmosphere white dwarfs
is changed (bottom row).

###  6.5 Velocity correction

As discussed in Sect. 2.2 and Sect. 3, by default our simulation includes the
age-velocity dispersion relation to compensate for a bias against high-
velocity objects in the 40 pc data set. To see if this correction has a
significant effect on the derived stellar formation history, we create a
simulation which does not include the age-velocity dispersion relation.
Instead, it assumes a constant scale height equal to 250 pc. This values does
not particularly matter because it only changes the absolute Galactic stellar
formation rate which is not directly constrained in this work. This new
simulation is then fitted to the observed 40 pc data. The best fit $\chi^{2}$
and $\tau_{\mathrm{max}}$ are shown in Table 2. These are in agreement with
the values found from the original simulation where the age-velocity
dispersion is included. The centre row of Fig. 11 shows the residuals and the
probability distribution of stellar formation history for this test. As can be
seen, there is no significant difference between the derived stellar formation
histories, indicating that the velocity dispersion correction does not affect
the derived stellar formation history significantly. This is in part because
the correlation between white dwarf total age and absolute $G$ magnitude is
weak, as the other important parameter, the progenitor mass, is left unchanged
in this test simulation.

###  6.6 Binary correction

Sect. 3.1 has described how the time delay correction arising from binary
mergers is included in our simulation. In this section we assess whether this
effect is significant for the derived stellar formation history. To do this, a
simulation is created where the binary time delay correction from Temmink et
al. 2020 is not included. As described in Sect. 3.1 in most cases the
inclusion of a binary correction means that more bright white dwarfs are
produced for the same stellar formation history. Thus, when fitting the 40 pc
sample, a simulation with no binary correction would have to produce less dim
white dwarfs which is done by decreasing $\tau_{\mathrm{max}}$ of the fitted
stellar formation history. The results are shown in Table 2 and Fig. 11. From
the $\chi^{2}$ information it is clear that the exclusion of merger time
delays does not change the best fit ECDF by more than the Poisson errors.

###  6.7 Fraction of He-atmosphere white dwarfs

In the original simulation, the fraction of He-atmosphere white dwarfs has
been set as 25% based on the 40 pc sample. To test whether this particular
choice has a significant effect on the derived stellar formation history, we
run two additional simulations where the fraction is set to 0% and to 50%. In
general, the effect of having more He-atmosphere white dwarfs is highly non-
linear, as for some white dwarf cooling ages and masses the $G$ magnitudes are
dimmer than H-atmosphere white dwarfs, whereas for others the $G$ magnitudes
are brighter.

In Table 2 and bottom row of Fig. 11 we show the results of this test. It is
clear that the particular choice of He-atmosphere white dwarf fraction has a
small effect on the ECDF or the derived SFH. However, we must note that the
observed ECDF is better fitted at dimmer magnitudes with a 0% He-fraction.
This is an interesting result as it suggests that having a well characterised
spectroscopic sample with precise atmospheric parameters and chemical
compositions may be as important as finding the right IMF and IFMR for fitting
the white dwarf absolute magnitude distribution.

###  6.8 White dwarf cooling tables

To test whether the white dwarf cooling tables have an effect on the derived
stellar formation history, we run a simulation with cooling tables generated
from the LPCODE stellar evolutionary code (Althaus et al., 2005; Althaus et
al., 2010b; Althaus et al., 2012; Rohrmann et al., 2012; Althaus et al., 2013;
Salaris et al., 2013; Camisassa et al., 2016; Miller Bertolami, 2016;
Camisassa et al., 2017; Rohrmann, 2018; Camisassa et al., 2019). In general,
the cooling tables directly affect the Gaia $G$ magnitude and the differences
between the Montréal tables and the LPCODE tables are complex. However, if one
table was to give consistently dimmer $G$ magnitudes, then the best fit values
of $\tau_{\mathrm{max}}$ would decrease when fitting the same sample of data.
For the same white dwarf mass and cooling ages, we find that the LPCODE
results in dimmer magnitudes.

Table 2 and Fig. 12 show the best fit values of $\chi^{2}$ and
$\tau_{\mathrm{max}}$ for this test. Based on $\chi^{2}$ this assumption
disagrees the most with the hypothesis of a uniform, one Galactic component
SFH. Therefore, a more complex SFH would be needed for the observed 40 pc
sample if the LPCODE is more accurate. However, it cannot be much different
from a uniform distribution, such as the SFH of Tremblay et al. (2014) or
Fantin et al. (2019), as the residuals for those are much larger (see Fig. 8).

Figure 12: This figure is similar to Figs. 10 and 11, but here we plot the
results of the tests where the white dwarf cooling tables are changed to those
generated by the LPCODE stellar evolutionary code.

###  6.9 Mass correction and double degenerate white dwarf removal

In this section, we switch from changing the different assumptions in the
simulation, and instead look at the observed 40 pc sample. In the previous
sections, we modified the original 40 pc white dwarf sample from McCleery et
al. (2020) and O’Brien et al. (2023) by applying a mass correction defined in
Table 1 for the low-mass problem of cool white dwarfs (see Sect. 2), and then
by removing any white dwarf whose mass was below the mass limit from single
star evolution, i.e. double degenerate candidates. A question arises whether
this modification of the 40 pc sample has an effect on the derived stellar
formation history. Therefore, to test this, we look at the following three
scenarios:

  1. 1.

The 40 pc sample used in all previous sections and described above, i.e. the
masses of white dwarfs have been corrected and double degenerates have been
removed. For this test we will refer to this sample as the "Corrected 40 pc
sample".

  2. 2.

The "Uncorrected 40 pc sample" refers to the 40 pc sample where mass
correction was not applied, but any white dwarf whose mass is below the limit
from single star evolution was removed.

  3. 3.

The last sample is the "Whole 40 pc sample" of McCleery et al. (2020) and
O’Brien et al. (2023), where no mass correction is applied and no white dwarf
has been removed.

In Fig. 13 we show the empirical cumulative distributions and the histograms
of the absolute $G$ magnitude for the three 40 pc samples. Over plotted on the
histogram are the Poissonian errors for the "Corrected 40 pc sample" as a
reference to judge whether the other two samples are significantly different.
Note that these errors are also plotted in Fig. 2.

It is clear from the top plot that the ECDFs of the "Uncorrected 40 pc sample"
is significantly different from the other two samples. As the positive mass
correction is not applied, the removal of white dwarfs with a mass below the
mass limit results in a larger number of white dwarfs removed when compared to
the "Corrected sample", i.e. 280 white dwarfs are removed as opposed to 110
which were removed to create the "Corrected 40 pc sample". We note that this
extra removal is likely spurious, i.e. not removing true double degenerates.
Nevertheless, these lower mass white dwarfs are cooler and older, and thus
their removal in the "Uncorrected sample" results in a histogram with more
bright white dwarfs. Predictably, this should mean a SFH which creates white
dwarfs closer to present time.

Figure 13: Top plot: A comparison of ECDFs for the 40 pc sample with different
mass cuts and corrections applied. The pink points represent the 40 pc sample
we have used in the rest of the paper, in this plot we refer to it as the
corrected 40 pc sample. This sample was corrected using the correction in
Table 1 and any white dwarf with a mass below the single star formation mass
limit was removed. The green triangles represent the sample where the mass
corrections were not applied, but any white dwarf with a mass below the single
star formation mass limit was removed. The black dots represent the complete
40 pc sample with no mass correction and no double degenerate removal, even
the white dwarfs whose mass is impossible if one assumes single star
formation. Bottom left plot: The histograms of the absolute $G$ magnitude for
the samples introduced in the top plot. The same colours are used to represent
the different samples as in the top plot. Bottom right plot: The best fitted
uniform, one-component stellar formation histories found from fitting the
different variations of the 40 pc sample. The same colours represent the
samples defined in top plot. This plot is similar to the left plots shown in
Fig. 10, 11 and 12.

In the bottom right plot of Fig. 13 we show the best fitted uniform, one-
component stellar formation histories for all three samples. In all cases the
simulation stays the same, only the observed sample is changed when fitting.
As before, we used MCMC with 1000 walkers to derive the probability density
functions of the best fitted SFHs. Table 2 shows the $\chi^{2}$ and
$\tau_{\mathrm{max}}$ of the fits as before. As expected, the "Corrected" and
the "Whole" samples result in very similar SFHs. However, the "Uncorrected
sample" results in a SFH which starts producing stars closer to present time,
since this sample has more bright white dwarfs when compared to the other two
samples. In all cases the reduced $\chi^{2}$ is close to 1, indicating that
uniform, one component SFH is a good fit to the data.

Overall, the "Corrected sample" is in better agreement with the population
synthesis expectations (i.e. white dwarfs cooling at constant mass) than the
"Uncorrected sample" (see Sect. 2 for more detailed discussion). Thus, in the
following we proceed with the stellar formation history derived from the
"Corrected 40 pc sample".

##  7 Discussion and conclusions

We have demonstrated that the white dwarf cooling tables and initial mass
function are the main assumptions that could allow for a deviation from
uniform stellar formation history with one galactic component. We conclude
that the white dwarf cooling models of Bédard et al. (2020) support a uniform
stellar formation history for the 40 pc sample, but this may not be the case
if other white dwarf cooling models such as the LPCODE were to be used.

By comparing our uniform, one Galactic component SFH with other SFHs we can
conclude that some minor $\approx$ 50% deviation from uniformity cannot be
ruled out, such as the SFH of Reid et al. (2007) or Mor et al. (2019).
However, more peaked SFHs of Tremblay et al. (2014) and Fantin et al. (2019)
are rejected by the data. The studies of Tremblay et al. (2014) and Fantin et
al. (2019) both use the same white dwarf cooling tables as us, albeit older
versions. As such, it is unlikely that the disagreement with those SFHs is due
to the cooling tables.

It is important to note that Sect. 4 has shown that increasing the number of
Galactic components could also improve the fit without needing to change other
astrophysical parameters, and still assuming separate, but uniform SFHs for
each component. Therefore, an argument can be made that any small deviation
from the constant formation history derived in this work can be compensated by
a small change in the input astrophysics. External constraints on these
parameters, such as accurate radial velocities for all stars in the sample,
better white dwarf cooling models or better constraints on the IMF, will be
needed to break this degeneracy.

Our conclusion is that the 40 pc sample of white dwarfs supports a uniform
stellar formation history given the current error bars on the various
astrophysical relations relevant for predicting the properties of the local
white dwarf population. Studies in need of a stellar formation history can
safely use the assumption of a constant formation history, with the onset of
stellar formation (age of the disc) and the number of Galactic components
within the local population, being the only parameters which remain difficult
to constrain.

The onset of stellar formation ($\tau_{\rm max}$) in our local region of the
Milky Way is linked to over two decades of work on the determination of the
age of the Galactic disc from white dwarf luminosity functions. Our main
result is that several significant uncertainties remain in the estimation of
this parameter. We have shown that changes in the IMF and white dwarf cooling
models can all change this parameter by up to $\pm 1$ Gyr. Furthermore, our
tests with different white dwarf cooling models is far from exhaustive.
Recently, it has been shown that our theoretical understanding of white dwarf
cooling is still lacking in some respects, as a missing cooling process
related to crystallisation has been identified when comparing theoretically
predicted and observed white dwarf luminosity functions (Tremblay et al.,
2019; Cheng et al., 2019; Cheng et al., 2020; Bauer et al., 2020; Blouin et
al., 2021; Camisassa et al., 2021). Therefore, we refrain from updating the
age of the Galactic disc based on our results, but note that our standard
simulation finds a value of 10.6 $\pm$ 0.5 Gyr, and that additional systematic
uncertainties may be as much as $\pm 1$ Gyr.

## Acknowledgements

The authors would like to thank the anonymous referee for their helpful
comments. EC and PET were supported by grant ST/T000406/1 from the Science and
Technology Facilities Council (STFC). PET has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant agreement No. 101020057). 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.

## Data Availability Statement

The observational data used in this article are published in McCleery et al.
(2020); Gentile Fusillo et al. (2021) and O’Brien et al. (2023). The results
derived in this article will be shared on a reasonable request to the
corresponding author.

## References

  * Althaus et al. (2005) Althaus L. G., García-Berro E., Isern J., Córsico A. H., 2005, A&A, 441, 689
  * Althaus et al. (2010a) Althaus L. G., Córsico A. H., Isern J., García-Berro E., 2010a, A&ARv, 18, 471
  * Althaus et al. (2010b) Althaus L. G., García-Berro E., Renedo I., Isern J., Córsico A. H., Rohrmann R. D., 2010b, ApJ, 719, 612
  * Althaus et al. (2012) Althaus L. G., García-Berro E., Isern J., Córsico A. H., Miller Bertolami M. M., 2012, A&A, 537, A33
  * Althaus et al. (2013) Althaus L. G., Miller Bertolami M. M., Córsico A. H., 2013, A&A, 557, A19
  * Andrews et al. (2015) Andrews J. J., Agüeros M. A., Gianninas A., Kilic M., Dhital S., Anderson S. F., 2015, ApJ, 815, 63
  * Anguiano et al. (2017) Anguiano B., Rebassa-Mansergas A., García-Berro E., Torres S., Freeman K. C., Zwitter T., 2017, MNRAS, 469, 2102
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Aumer & Binney (2009) Aumer M., Binney J. J., 2009, MNRAS, 397, 1286
  * Bailer-Jones et al. (2018) Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R., 2018, AJ, 156, 58
  * Bailer-Jones et al. (2021) Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Demleitner M., Andrae R., 2021, AJ, 161, 147
  * Bauer et al. (2020) Bauer E. B., Schwab J., Bildsten L., Cheng S., 2020, ApJ, 902, 93
  * Bédard et al. (2020) Bédard A., Bergeron P., Brassard P., Fontaine G., 2020, ApJ, 901, 93
  * Belokurov (2013) Belokurov V., 2013, New Astron. Rev., 57, 100
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Bensby et al. (2014) Bensby T., Feltzing S., Oey M. S., 2014, A&A, 562, A71
  * Bergeron et al. (2019) Bergeron P., Dufour P., Fontaine G., Coutu S., Blouin S., Genest-Beaulieu C., Bédard A., Rolland B., 2019, ApJ, 876, 67
  * Besla et al. (2010) Besla G., Kallivayalil N., Hernquist L., van der Marel R. P., Cox T. J., Kereš D., 2010, ApJ, 721, L97
  * Bessell & Stringfellow (1993) Bessell M. S., Stringfellow G. S., 1993, ARA&A, 31, 433
  * Blouin et al. (2021) Blouin S., Daligault J., Saumon D., 2021, ApJ, 911, L5
  * Bond et al. (2010) Bond N. A., et al., 2010, ApJ, 716, 1
  * Bovy et al. (2012a) Bovy J., Rix H.-W., Hogg D. W., 2012a, ApJ, 751, 131
  * Bovy et al. (2012b) Bovy J., Rix H.-W., Liu C., Hogg D. W., Beers T. C., Lee Y. S., 2012b, ApJ, 753, 148
  * Buckner & Froebrich (2014) Buckner A. S. M., Froebrich D., 2014, MNRAS, 444, 290
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Camisassa et al. (2016) Camisassa M. E., Althaus L. G., Córsico A. H., Vinyoles N., Serenelli A. M., Isern J., Miller Bertolami M. M., García–Berro E., 2016, ApJ, 823, 158
  * Camisassa et al. (2017) Camisassa M. E., Althaus L. G., Rohrmann R. D., García-Berro E., Torres S., Córsico A. H., Wachlin F. C., 2017, ApJ, 839, 11
  * Camisassa et al. (2019) Camisassa M. E., et al., 2019, A&A, 625, A87
  * Camisassa et al. (2021) Camisassa M. E., Althaus L. G., Torres S., Córsico A. H., Rebassa-Mansergas A., Tremblay P.-E., Cheng S., Raddi R., 2021, A&A, 649, L7
  * Casagrande et al. (2011) Casagrande L., Schönrich R., Asplund M., Cassisi S., Ramírez I., Meléndez J., Bensby T., Feltzing S., 2011, A&A, 530, A138
  * Cheng et al. (2019) Cheng S., Cummings J. D., Ménard B., 2019, ApJ, 886, 100
  * Cheng et al. (2020) Cheng S., Cummings J. D., Ménard B., Toonen S., 2020, ApJ, 891, 160
  * Cole et al. (2000) Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168
  * Cukanovaite et al. (2021) Cukanovaite E., Tremblay P.-E., Bergeron P., Freytag B., Ludwig H.-G., Steffen M., 2021, MNRAS, 501, 5274
  * Cummings et al. (2016) Cummings J. D., Kalirai J. S., Tremblay P. E., Ramirez-Ruiz E., 2016, ApJ, 818, 84
  * Cummings et al. (2018) Cummings J. D., Kalirai J. S., Tremblay P. E., Ramirez-Ruiz E., Choi J., 2018, ApJ, 866, 21
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * De Lucia et al. (2006) De Lucia G., Springel V., White S. D. M., Croton D., Kauffmann G., 2006, MNRAS, 366, 499
  * Du et al. (2018) Du C., Li H., Liu S., Donlon T., Newberg H. J., 2018, ApJ, 863, 87
  * El-Badry et al. (2018) El-Badry K., Rix H.-W., Weisz D. R., 2018, ApJ, 860, L17
  * Fantin et al. (2019) Fantin N. J., et al., 2019, ApJ, 887, 148
  * Fantin et al. (2021) Fantin N. J., et al., 2021, ApJ, 913, 30
  * Feltzing & Gonzalez (2001) Feltzing S., Gonzalez G., 2001, A&A, 367, 253
  * Fontaine et al. (2001) Fontaine G., Brassard P., Bergeron P., 2001, PASP, 113, 409
  * Frebel & Norris (2015) Frebel A., Norris J. E., 2015, ARA&A, 53, 631
  * Gaia Collaboration et al. (2016a) Gaia Collaboration et al., 2016a, A&A, 595, A1
  * Gaia Collaboration et al. (2016b) Gaia Collaboration et al., 2016b, A&A, 595, A2
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Gallart et al. (2019) Gallart C., Bernard E. J., Brook C. B., Ruiz-Lara T., Cassisi S., Hill V., Monelli M., 2019, Nature Astronomy, 3, 932
  * García-Berro & Oswalt (2016) García-Berro E., Oswalt T. D., 2016, New Astron. Rev., 72, 1
  * Garcia-Berro et al. (1988) Garcia-Berro E., Hernanz M., Isern J., Mochkovitch R., 1988, Nature, 333, 642
  * García-Berro et al. (2010) García-Berro E., et al., 2010, Nature, 465, 194
  * Gentile Fusillo et al. (2019) Gentile Fusillo N. P., et al., 2019, MNRAS, 482, 4570
  * Gentile Fusillo et al. (2021) Gentile Fusillo N. P., et al., 2021, MNRAS, 508, 3877
  * Gilmore et al. (1989) Gilmore G., Wyse R. F. G., Kuijken K., 1989, ARA&A, 27, 555
  * Grillmair & Carlin (2016) Grillmair C. J., Carlin J. L., 2016, in Newberg H. J., Carlin J. L., eds, Astrophysics and Space Science Library Vol. 420, Tidal Streams in the Local Group and Beyond. p. 87 (arXiv:1603.08936), doi:10.1007/978-3-319-19336-6_4
  * Guedes et al. (2013) Guedes J., Mayer L., Carollo M., Madau P., 2013, ApJ, 772, 36
  * Harris et al. (2006) Harris H. C., et al., 2006, AJ, 131, 571
  * Hayden et al. (2015) Hayden M. R., et al., 2015, ApJ, 808, 132
  * Hayden et al. (2017) Hayden M. R., Recio-Blanco A., de Laverny P., Mikolaitis S., Worley C. C., 2017, A&A, 608, L1
  * Haywood et al. (2013) Haywood M., Di Matteo P., Lehnert M. D., Katz D., Gómez A., 2013, A&A, 560, A109
  * Haywood et al. (2018) Haywood M., Di Matteo P., Lehnert M. D., Snaith O., Khoperskov S., Gómez A., 2018, ApJ, 863, 113
  * Heintz et al. (2022) Heintz T. M., Hermes J. J., El-Badry K., Walsh C., van Saders J. L., Fields C. E., Koester D., 2022, ApJ, 934, 148
  * Helmi (2008) Helmi A., 2008, A&ARv, 15, 145
  * Helmi (2020) Helmi A., 2020, ARA&A, 58, 205
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Hollands et al. (2018) Hollands M. A., Tremblay P. E., Gänsicke B. T., Gentile-Fusillo N. P., Toonen S., 2018, MNRAS, 480, 3942
  * Hou & Han (2014) Hou L. G., Han J. L., 2014, A&A, 569, A125
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Ibata et al. (2002) Ibata R. A., Lewis G. F., Irwin M. J., Quinn T., 2002, MNRAS, 332, 915
  * Isern (2019) Isern J., 2019, ApJ, 878, L11
  * Isern et al. (1998) Isern J., García-Berro E., Hernanz M., Mochkovitch R., Torres S., 1998, ApJ, 503, 239
  * Jimenez et al. (1998) Jimenez R., Flynn C., Kotoneva E., 1998, MNRAS, 299, 515
  * Jofré & Weiss (2011) Jofré P., Weiss A., 2011, A&A, 533, A59
  * Kawata & Chiappini (2016) Kawata D., Chiappini C., 2016, Astronomische Nachrichten, 337, 976
  * Kereš et al. (2005) Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2
  * Kilic et al. (2017) Kilic M., Munn J. A., Harris H. C., von Hippel T., Liebert J. W., Williams K. A., Jeffery E., DeGennaro S., 2017, ApJ, 837, 162
  * Koester & Chanmugam (1990) Koester D., Chanmugam G., 1990, Reports on Progress in Physics, 53, 837
  * Koppelman et al. (2018) Koppelman H., Helmi A., Veljanoski J., 2018, ApJ, 860, L11
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kroupa (2002) Kroupa P., 2002, Science, 295, 82
  * Kruijssen et al. (2019) Kruijssen J. M. D., Pfeffer J. L., Reina-Campos M., Crain R. A., Bastian N., 2019, MNRAS, 486, 3180
  * Lindegren et al. (2021a) Lindegren L., et al., 2021a, A&A, 649, A2
  * Lindegren et al. (2021b) Lindegren L., et al., 2021b, A&A, 649, A4
  * Liu & Chaboyer (2000) Liu W. M., Chaboyer B., 2000, ApJ, 544, 818
  * López-Sanjuan et al. (2022) López-Sanjuan C., et al., 2022, A&A, 658, A79
  * Majewski (1992) Majewski S. R., 1992, ApJS, 78, 87
  * McCleery et al. (2020) McCleery J., et al., 2020, MNRAS, 499, 1890
  * Miller Bertolami (2016) Miller Bertolami M. M., 2016, A&A, 588, A25
  * Minchev & Famaey (2010) Minchev I., Famaey B., 2010, ApJ, 722, 112
  * Minchev et al. (2015) Minchev I., Martig M., Streich D., Scannapieco C., de Jong R. S., Steinmetz M., 2015, ApJ, 804, L9
  * Mo et al. (1998) Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319
  * Mochkovitch et al. (1990) Mochkovitch R., Garcia-Berro E., Hernanz M., Isern J., Panis J. F., 1990, A&A, 233, 456
  * Mor et al. (2019) Mor R., Robin A. C., Figueras F., Roca-Fàbrega S., Luri X., 2019, A&A, 624, L1
  * Myeong et al. (2019) Myeong G. C., Vasiliev E., Iorio G., Evans N. W., Belokurov V., 2019, MNRAS, 488, 1235
  * Napiwotzki et al. (2020) Napiwotzki R., et al., 2020, A&A, 638, A131
  * Niederste-Ostholt et al. (2010) Niederste-Ostholt M., Belokurov V., Evans N. W., Peñarrubia J., 2010, ApJ, 712, 516
  * Nissen & Schuster (2010) Nissen P. E., Schuster W. J., 2010, A&A, 511, L10
  * Noh & Scalo (1990) Noh H.-R., Scalo J., 1990, ApJ, 352, 605
  * Norris & Ryan (1991) Norris J. E., Ryan S. G., 1991, ApJ, 380, 403
  * O’Brien et al. (2023) O’Brien M. W., et al., 2023, MNRAS, 518, 3055
  * Odenkirchen et al. (2001) Odenkirchen M., et al., 2001, ApJ, 548, L165
  * Pauli et al. (2006) Pauli E. M., Napiwotzki R., Heber U., Altmann M., Odenkirchen M., 2006, A&A, 447, 173
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A20
  * Portegies Zwart & Verbunt (1996) Portegies Zwart S. F., Verbunt F., 1996, A&A, 309, 179
  * Raddi et al. (2022) Raddi R., et al., 2022, A&A, 658, A22
  * Rebassa-Mansergas et al. (2015) Rebassa-Mansergas A., et al., 2015, MNRAS, 450, 743
  * Rebassa-Mansergas et al. (2021) Rebassa-Mansergas A., et al., 2021, MNRAS, 505, 3165
  * Reid et al. (2007) Reid I. N., Turner E. L., Turnbull M. C., Mountain M., Valenti J. A., 2007, ApJ, 665, 767
  * Rohrmann (2018) Rohrmann R. D., 2018, MNRAS, 473, 457
  * Rohrmann et al. (2012) Rohrmann R. D., Althaus L. G., García-Berro E., Córsico A. H., Miller Bertolami M. M., 2012, A&A, 546, A119
  * Rowell (2013) Rowell N., 2013, MNRAS, 434, 1549
  * Salaris & Weiss (2002) Salaris M., Weiss A., 2002, A&A, 388, 492
  * Salaris et al. (2013) Salaris M., Althaus L. G., García-Berro E., 2013, A&A, 555, A96
  * Sellwood & Binney (2002) Sellwood J. A., Binney J. J., 2002, MNRAS, 336, 785
  * Siebert et al. (2008) Siebert A., et al., 2008, MNRAS, 391, 793
  * Siegert (2019) Siegert T., 2019, A&A, 632, L1
  * Smith et al. (2009) Smith M. C., et al., 2009, MNRAS, 399, 1223
  * Soderblom (2010) Soderblom D. R., 2010, ARA&A, 48, 581
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Speagle et al. (2014) Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15
  * Temmink et al. (2020) Temmink K. D., Toonen S., Zapartas E., Justham S., Gänsicke B. T., 2020, A&A, 636, A31
  * Tononi et al. (2019) Tononi J., Torres S., García-Berro E., Camisassa M. E., Althaus L. G., Rebassa-Mansergas A., 2019, A&A, 628, A52
  * Toonen et al. (2012) Toonen S., Nelemans G., Portegies Zwart S., 2012, A&A, 546, A70
  * Torres et al. (2019) Torres S., Cantero C., Rebassa-Mansergas A., Skorobogatov G., Jiménez-Esteban F. M., Solano E., 2019, MNRAS, 485, 5573
  * Torres et al. (2021) Torres S., Rebassa-Mansergas A., Camisassa M. E., Raddi R., 2021, MNRAS, 502, 1753
  * Tremblay et al. (2011) Tremblay P. E., Bergeron P., Gianninas A., 2011, ApJ, 730, 128
  * Tremblay et al. (2014) Tremblay P. E., Kalirai J. S., Soderblom D. R., Cignoni M., Cummings J., 2014, ApJ, 791, 92
  * Tremblay et al. (2016) Tremblay P. E., Cummings J., Kalirai J. S., Gänsicke B. T., Gentile-Fusillo N., Raddi R., 2016, MNRAS, 461, 2100
  * Tremblay et al. (2019) Tremblay P.-E., et al., 2019, Nature, 565, 202
  * Tremblay et al. (2020) Tremblay P. E., et al., 2020, MNRAS, 497, 130
  * Van Rossum & Drake (2009) Van Rossum G., Drake F. L., 2009, Python 3 Reference Manual. CreateSpace, Scotts Valley, CA
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wegg & Phinney (2012) Wegg C., Phinney E. S., 2012, MNRAS, 426, 427
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Winget et al. (1987) Winget D. E., Hansen C. J., Liebert J., van Horn H. M., Fontaine G., Nather R. E., Kepler S. O., Lamb D. Q., 1987, ApJ, 315, L77
  * del Peloso et al. (2005a) del Peloso E. F., da Silva L., Porto de Mello G. F., 2005a, A&A, 434, 275
  * del Peloso et al. (2005b) del Peloso E. F., da Silva L., Porto de Mello G. F., Arany-Prado L. I., 2005b, A&A, 440, 1153