11institutetext: Department of Astrophysics, University of Vienna,
Türkenschanzstrasse 17, 1180 Wien, Austria
11email<EMAIL_ADDRESS>
# An empirical isochrone archive for nearby open clusters

Alena Rottensteiner 11 Stefan Meingast 11

(Received 10 August 2023 / Accepted 3 June 2024)

###### Abstract

Context. The ages of star clusters and co-moving stellar groups contain
essential information about the Milky Way. Their special properties and
placement throughout the galactic disk make them excellent tracers of galactic
structure and key components to unlocking its star formation history. Yet,
even though the importance of stellar population ages has been widely
recognized, their determination remains a challenging task often associated
with highly model-dependent and uncertain results.

Aims. We propose a new approach to this long-standing problem, which relies on
empirical isochrones of known clusters extracted from high-quality
observational data. These purely observation-based data products open up the
possibility of relative age determination, free of stellar evolution model
assumptions.

Methods. For the derivation of the empirical isochrones, we used a combination
of the statistical analysis tool principal component analysis for
preprocessing and the supervised machine learning method support vector
regression for curve extraction. To improve the statistical reliability of our
result, we defined the empirical isochrone of a color-magnitude diagram (CMD)
of a cluster as the median calculated from a set of $n_{\mathrm{boot}}=1000$
curves derived from bootstrapped data. The algorithm requires no physical
priors, is computationally fast, and can easily be generalized over a large
range of CMD combinations and evolutionary stages of clusters.

Results. We provide empirical isochrones in all _Gaia_ DR2 and DR3 color
combinations for 83 nearby clusters ($d<500$ pc), which cover an estimated age
range of 7 Myr to 3 Gyr. In doing so, we pave the way for a relative
comparison between individual stellar populations based on an age-scaling
ladder of empirical isochrones of known clusters. Furthermore, due to the
exceptional precision of the available observational data, we report accurate
lower main sequence empirical isochrones for many clusters in our sample,
which are of special interest as this region is known to be especially complex
to model. We validate our method and results by comparing the extracted
empirical isochrones to cluster ages in the literature. We also investigate
the added information that empirical isochrones covering the lower main
sequence can provide on case studies of the IC 4665 cluster and the Meingast 1
stream.

Conclusions. The archive of empirical isochrones offers a novel approach to
validating age estimates and can be used as an age-scaling ladder or age
brackets for new populations and serve as calibration data for further
constraining stellar evolution models.

###### Key Words.:

Stars: Hertzsprung-Russel and C-M diagrams - Methods: data analysis - Methods:
statistical - Open clusters and associations: general

##  1 Introduction

Among the wide variety of sources that can be observed in the plane of the
night sky, star clusters have long been held in a particularly prominent
regard. This fascination is at least partly due to the fact that these complex
systems, comprising tens, hundreds, or even thousands of stars that formed in
the same molecular cloud, often map spectacular constellations onto the
firmament. But also besides their picturesque appearances, star clusters have
long been established as objects of dedicated study in various fields, among
them the research into the formation and evolution of galactic structure
across all scales. They are exceptional laboratories for studying
characteristics of a stellar sample with largely homogeneous initial
conditions and formation history, and especially for open clusters, their
predominant spatial occurrence within the galactic plane identifies them as
strategic targets for probing different galactic properties. Their generally
young ages and formation history forge a deep connection between clusters and
past or still active star formation sites.

The advent of extensive space-based, all-sky astrometric surveys, such as
Hipparcos (Perryman et al. 1997) and later the _Gaia_ mission with its various
data releases (Gaia Collaboration et al. 2016, 2017, 2018b, 2022) and
unprecedented observational accuracy and quantity, has heralded new insights
not only into the shape and structure of clusters themselves (e.g., Meingast
et al. 2021, cluster coronae) but also regarding their mean properties (e.g.,
Bossini et al. 2019; Cantat-Gaudin et al. 2020; Dias et al. 2021). Moreover,
access to parallaxes, proper motions, and radial velocities has caused the
scientific community to move away from studying only the clusters visible in
positional space and instead focus on exploring the combined spatial and
kinematic phase space. As a result, an ever-growing number of clusters and
stellar groups is being discovered using various methods measuring
overdensities within this space (e.g., Cantat-Gaudin et al. 2018; Cantat-
Gaudin & Anders 2020; Kerr et al. 2021, 2022a, 2022b; Castro-Ginard et al.
2022; Ratzenböck et al. 2022; Hunt & Reffert 2023), while at the same time
conglomerates formerly regarded as clusters have been revealed to be only
asterisms (e.g., Cantat-Gaudin et al. 2018; Cantat-Gaudin & Anders 2020). Both
findings are crucial for advancing the cluster and association census in the
Milky Way. This, in turn, has consequences for our understanding of spatial
structures on a galactic scale, as shown, for example, by Castro-Ginard et al.
(2021), who successfully traced the spiral arm structure of the Milky Way
using open clusters.

Open clusters also greatly aid the accessibility of our galactic past due to
their close relation to star-forming regions and as direct representations of
star formation history. However, to fully access the temporal axis of
structure formation and evolution on larger-scale objects, one requires
knowledge about a crucial yet elusive physical parameter – the age of a
cluster or a stellar population. The age parameter plays a critical role in
developing timescales, for instance, for formation processes of large-scale
star formation regions. For example, Ratzenböck et al. (2023) recently mapped
a possible evolutionary scenario for the nearest active star-forming site,
Scorpius-Centaurus (Sco-Cen), using only the information deduced from the over
thirty different star clusters they isolated in the 6D phase space. Cluster
ages are also a vital piece for enhancing our understanding of molecular cloud
formation and dispersal (e.g., Murray 2011; Adamo et al. 2013; Chevance et al.
2020), or even for constraining possible windows for exoplanet formation
(David & Hillenbrand 2015).

However, ages are not among the list of directly accessible parameters of
stars, which means that, despite their importance across many different
physical scales, their determination remains a notoriously difficult task to
this day. This is especially true for stellar population ages. Various age
estimation techniques have been developed for these objects, but their results
are seldom unanimous. Even for well-studied, local clusters, as for example
the Pleiades (Melotte 22), for which high-quality observational data are
available, the age estimation process is often characterized by the need for
prior assumptions and numerical models, leading to, at times, vastly deviating
age estimates (see Fig. 1). The heavy model dependence of population age
estimation also holds for empirical methods to some extent, as they need to be
calibrated using benchmark clusters and, therefore, partially inherit their
associated uncertainties. Notably, the Pleiades cluster seen in Fig. 1 with
its evident discrepancies in age estimates, is among the most commonly used
benchmarks (e.g., Bell et al. 2012; Olivares et al. 2018a; Messina et al.
2022; Brandner et al. 2023). Moreover, uncertainties or gaps in current
stellar evolution theories, computational limitations in the implementation of
stellar interior physics, as well as measurement errors in observational data
cause further complications and together preclude the determination of
absolute ages with universal, unbiased, and academically undisputed methods
(see, e.g., Soderblom 2010, and references therein).

To circumvent many of the drawbacks associated with the age determination of
stellar populations, in particular regarding stellar evolutionary model
dependencies, we have developed an innovative approach to age estimation:
Using only the information from observational Hertzsprung-Russel diagrams
(HRDs) of stellar populations, we extract empirical isochrones for open
clusters and thus pave the way for the derivation of relative cluster ages.
Their precise representation of the source distribution can add vital
information about the lower mass end of a stellar population and could be
valuable for better constraining the physical models of stellar interiors in
the low-mass regime. Based exclusively on widely available observational data,
the distinguishing property of our approach, compared to other age estimation
methods, is its purely factual nature. By providing empirical isochrones for a
comprehensive, representative, and relevant selection of nearby open clusters,
we establish a way of quantifying ages via a relative and homogeneous
comparison between individual stellar populations.

The remainder of this work is structured as follows: We start by briefly
explaining the current state of the art of stellar population age
determination and the complications arising from different methods, that
inspired the creation of the empirical isochrone archive in Sect. 2. We
describe the data and cluster memberships we used in Sect. 3 before reviewing
our methods and the isochrone extraction algorithm in Sect. 4. The results are
presented in Section 5 and validated against age estimates found in
contemporary literature. We follow up with a discussion of the isochrone
quality and significance, especially for the region of interest in the lower
main sequence, with the aid of case studies in Section 6. Lastly, we summarize
our findings in Sect. 7.

##  2 The current state of population age estimation

Since there is no measurable parameter that enables unequivocal age
determination, various methods have been developed to estimate stellar
population ages using different indicators. These include asteroseismology
(e.g., Pamos Ortega et al. 2022), evaporation ages (Pelkonen et al. 2024),
gyrochronology (e.g., Barnes 2007; Curtis et al. 2019), kinematic and
traceback ages (e.g., Miret-Roig et al. 2020; Kerr et al. 2022a), lithium
depletion (boundary) measurements (e.g., Basri et al. 1996; Binks & Jeffries
2014; Martín et al. 2018), or lithium equivalent widths (e.g., Jeffries et al.
2023). However, these methods all share a “specialist” limitation, in that
they may only be suited for specific age, temperature, or mass ranges, require
precise, work-intensive observed parameters, or can only be applied to either
individual stars or ensembles. (see, e.g., Table 1 in Soderblom 2010).

Contrasting these rather specialized methods with narrow application ranges is
an approach often praised for its generalistic nature – isochrone fitting.
Based on variants of the Hertzsprung-Russel diagram, particularly the color-
absolute magnitude diagram (CMD), it is a popular estimation technique with
only minimal observational data requirements. As a result of _Gaia_ providing
photometry and astrometry for billions of sources, isochrone fitting has
become perhaps the most used approach to age determination to date. Because of
its connection to the HRD, the method is widely applicable across stellar age
and mass ranges, for example from 0.5 Myr to 10 Gyr for
$M_{*}=0.01-1.4\leavevmode\nobreak\ M_{\odot}$ (Baraffe et al. 2015) or 10 Myr
to 12.5 Gyr for $M_{*}=0.09-14\leavevmode\nobreak\ M_{\odot}$ (Nguyen et al.
2022). There are two major drawbacks to fitting isochrones, though: On the one
hand, it heavily relies on the specific physics considered when modeling the
theoretical isochrones. As the different physical processes of stellar
evolution are highly complex and vary between different stellar ages and
masses (de Boer & Seggewiss 2008), this translates into a large number of free
parameters, which can produce systematic errors in age estimates across
different stellar evolution models and prior assumptions. Given the many
different available stellar evolution codes for model isochrones, such as
BHAC15 (Baraffe et al. 2015), DSEP-magnetic (Feiden 2016), BT-Settl (Allard
2014), PARSEC (Bressan et al. 2012; Marigo et al. 2017; Nguyen et al. 2022),
MIST (Choi et al. 2016; Dotter 2016) and BaSTI (Hidalgo et al. 2018;
Pietrinferni et al. 2021; Salaris et al. 2022), this dependence has a profound
impact on age estimates found in literature. On the other hand, the fitting
technique of the isochrones to the CMD data itself can strongly impact the
results. From the observational side, measurement errors and incomplete or
contaminated cluster member data are drivers of age uncertainties (de Boer &
Seggewiss 2008).

Additionally, observational cluster CMDs include an unknown amount of
unresolved binary stars. The extent of this contaminating population has been
estimated to range from about 29 to 50 % (Li & Shao 2022), and equal mass
unresolved binaries reside at the maximum possible displacement of
$-0.753\text{\,}\mathrm{m}\mathrm{a}\mathrm{g}$ above the main sequence.
Additionally, reddening effects of non-equal mass unresolved binaries have
been found to cause a shift of the corresponding data points closer, or right
onto, the equal mass binary sequence (e.g., Hurley & Tout 1998; Li et al.
2020). Consequently, it often appears disproportionately enlarged compared to
the expected number of equal mass binaries. The presence of such a well-
defined locus of contaminants harbors another challenge to isochronal age
determination. However, the impact of the binary sequence on the empirical
isochrones presented in this work has been specifically addressed and
quantified (see Appendix E.

Consequently, it is challenging to verify age estimates obtained through
different age estimation methods or to establish a uniform age for an object
using the literature (see, e.g., the comparison tables in Mamajek & Bell 2014;
Miret-Roig et al. 2020; Kerr et al. 2022b). In the following paragraphs, we
examine this situation’s implications, focusing on isochrone fitting. For a
more detailed explanation of the age determination methods discussed here,
readers can refer to the review works of Soderblom (2010); Soderblom et al.
(2014).

###  2.1 Age ambiguity

The apparent ambiguity of stellar ages derived with different methods is at
least in part caused by their use of different age indicators. For instance,
dynamical traceback ages, which measure the time since a group of stars was
most concentrated in space, seem to systematically underestimate the ages
compared to other techniques, such as isochrone fitting or lithium depletion
boundary measurements, of young stellar agglomerates ($\lesssim 50$ Myr). A
possible explanation for this observation could be that they already existed
for a few Myr without significant expansion, which would not be reflected in
the kinematic age indicator (Miret-Roig et al. 2024). Another example concerns
very young clusters, whose lithium depletion boundary ages have been reported
to systematically be about 50 % greater than those determined from their main
sequence turn-off points (Soderblom 2010, and references therein).

In contrast to methods based on different age indicators, the strong
dependence of the results on the individual model assumptions causes isochrone
fitting to produce significant age differences even though the overarching
determination technique remains the same. This seriously impacts the
reliability and the reproducibility of isochronal ages, as highlighted by Kerr
et al. (2022a, b). They compare isochronal ages from three stellar evolution
codes, namely PARSEC (Chen et al. 2015), BHAC15 (Baraffe et al. 2015), and
DSEP-Magnetic (Feiden 2016), on different stellar groups and find age
disparities of up to a factor of $\sim 2.6$ between the different models.

Figure 1: Comparison between two different isochronal ages for the Pleiades
cluster. The data points correspond to the _Gaia_ DR2 cluster selection of
Cantat-Gaudin & Anders (2020) , and the isochrones were produced with PARSEC,
using the parameters calculated in the respective works (Bossini et al. 2019;
Dias et al. 2021), indicated by B19 and D21 in the figure legend,
respectively.

Further investigations show that even employing theoretical isochrones from
the same model and the same member selection can lead to extreme age
differences, as in the case of the Pleiades cluster depicted in Fig. 1. The
plot shows the DR2 member selection of the cluster by Cantat-Gaudin & Anders
(2020)111In the original selection of the authors, an artifact appears in the
form of a vertical line of data points at the color index zero. It disappears
when cross-matching their membership list with DR3 data., along with
theoretical isochrones based on two published sets of age, metallicity, and
extinction values, which were determined via interpolation of a parameter grid
of PARSEC V1.2S isochrones. Both the theoretical isochrones displayed in the
figure (Bossini et al. 2019; Dias et al. 2021, denoted as B19 and D21)222The
authors used different parameter grids and _Gaia_ passband corrections. They
report a similar cluster extinction ($A_{V,\leavevmode\nobreak\
\textrm{Dias}}=0.168$, $A_{V,\leavevmode\nobreak\ \textrm{Bossini}}=0.14$);
however, Dias et al. (2021) treated the cluster metallicity as a free
parameter (Z${}_{\textrm{Dias}}$ = 0.032), while Bossini et al. (2019) assumed
solar metallicity for the Pleiades. present a good fit to the observations
over almost the entire populated dynamical range between ca. 0 and 10 absolute
G-band magnitudes. However, their absolute age difference amounts to 45 Myr,
which is around a 50 % deviation of the lower age estimate. Yet, they are
almost indistinguishable from each other, with the exception of the faint mass
end of the sequence ($\mathrm{M}_{\mathrm{G}}\gtrsim 11$ mag), where neither
isochrone presents a good fit to the data, and the the upper main sequence,
which lacks observations.

The fact that such inconclusive age estimates exist even for one of the
nearest, best-studied open clusters in the solar neighborhood, which is
moreover a popular benchmark for empirical estimation methods, illustrates the
extent of uncertainty that can be accrued when using isochrone fitting.

###  2.2 Isochrone blindspot

Figure 2: Depiction of the largely featureless region, dubbed the “isochrone
blindspot.” It covers an age span of approximately $100-500$ Myr and roughly
$0-12$ G magnitudes in the _Gaia_ passbands. The plot features PARSEC
isochrones (Bressan et al. 2012) with default parameters. The _dashed line_
indicates the upper main sequence region where observed sources are rare,
while the _dash-dotted line_ delineates the lower main sequence, which is a
region of interest regarding the information gained through empirical
isochrones.

As shown in the example of the Pleiades cluster, theoretical isochrones can
appear almost indiscernible over a large portion of the main sequence, despite
being associated with significantly deviating ages (see also, e.g., de Boer &
Seggewiss 2008; Squicciarini & Bonavita 2022). This effect is particularly
evident for clusters in an age range around $\sim 100-500$ Myr and is
prevalent across several photometric bands, not only for the specific set of
_Gaia_ passbands. Moreover, it is not even restricted to the optical
wavelength range, but can for example also be observed in the infrared $J$,
$H$ and $K_{S}$ filters. It should be noted, that the dynamical range affected
by the blindspot varies depending on the passbands used in the CMD.

Regarding the _Gaia_ photometric system, the overlapping of the theoretical
isochrones commences as soon as their main sequences333In this context we
refer to the lower main sequence as roughly $\mathrm{M}_{\mathrm{G}}\gtrsim
10$ mag, and to the upper main sequence when $\mathrm{M}_{\mathrm{G}}\lesssim
2$ mag. do no longer appear flatter and shifted to redder colors like those of
very young clusters ($\lesssim 30$ Myr). To illustrate this issue, Fig. 2
shows a set of PARSEC isochrones calculated with default metallicity and
extinction values for this age range in the three _Gaia_ DR3 color variations.
The ambiguity of the main sequence across large parts of the CMD requires the
isochronal age determination to hinge on available sources located in the
upper main sequence and beyond (_dashed line_), where the isochrones start to
deviate again. However, the main sequence turn-off and evolved star regions
only begin to get well populated for very evolved clusters, meaning no
reliable information for the depicted age range is usually available for
isochrone fitting. In addition, according to simulations and observational
trends, massive stars are more likely part of binary or multiple systems
(e.g., De Rosa et al. 2014, and references therein). If such a system was
observationally unresolved, the age determination could be seriously
influenced. The same holds for possible field star contaminants around the
turn-off region, simply because of the typically low number of observations in
this crucial region for isochrone fitting.

In short, from the point at which young clusters are no longer separated from
older populations via their shifted, flatter main sequence until the time at
which the main sequence turn-off has shifted to lower stellar masses and
becomes more populated for evolved clusters, there are no decisive features in
the shapes of theoretical isochrones over a sizeable dynamical range –
resulting in the appearance of an “isochrone blindspot.” Physically, this
blindspot is a consequence of the very stable and long-lasting burning
processes occurring during the main sequence phase of lower mass stars (e.g.
Barnes 2007; Squicciarini & Bonavita 2022), that results in only minute
variations of the luminosity and temperature, or color, over their main
sequence lifetimes. What is more, the blindspot does not only exist for
theoretical isochrones but also partly occurs when considering other age
determination methods for stellar populations. For instance, gyrochronology
only works for single, solar-type stars, whereas lithium depletion boundary
measurements or kinematic modeling apply only to young groups of $\lesssim
200$ Myr and $\lesssim 30$ Myr, respectively (Soderblom et al. 2014).

A possible workaround of the blindspot problem would be to consult the lower
main sequences of clusters, below the _dash-dotted line_ in Fig. 2. As can be
seen, for this region the model isochrones are again better distinguishable;
however, as depicted in the lower right panel of Fig. 1 for the Pleiades,
models often do not capture the observations in enough detail to permit an age
estimation using only this part of the CMD with theoretical isochrone fitting.
Empirical isochrones, on the other hand, very precisely match the shape of the
lower main sequence of a CMD, provided there are observations of faint stars
available for a given cluster. Comparing empirical lower main sequence
isochrones thus provides an additional piece of information to be leveraged
when deciding on an absolute age estimate, or can even be sufficient for a
relative bracketing or ordering of observed stellar population ages (Sect.
6.2).

###  2.3 Empirical isochrones – A largely untapped potential

Despite current limitations in technology and physical knowledge concerning
stellar evolution modeling, the recent rise in precision and quantity of
stellar photometry and astrometry is opening up a new, data-driven approach
for using isochrones as a decisive age indication tool: in an empirical form.
Both empirical and theoretical isochrones have a fundamental connection to the
HRD, which renders almost all aspects of stellar evolution accessible from a
2D representation. Consequently, empirical isochrones do not suffer from the
problem of being too “specialist” compared to other age determination methods
discussed in this section. However, in contrast to theoretical isochrones, the
empirical curves directly originate from the CMD, meaning they are purely
observation-based. Thus, in creating the empirical isochrone archive, we avoid
the caveats and drawbacks associated explicitly with the model dependency of
theoretical isochrones, while still taking advantage of the extensive
applicability of the constructs of isochrones in general. Since the CMD
comprises only three easily accessible observational parameters, empirical
isochrone extraction profits from a tiny list of requirements and no physical
complexity compared to the models necessary for computing theoretical
isochrones.

The value of empirical isochrones has already been recognized among the
scientific community, and (semi)-empirical isochrones have been extracted
using different techniques and published for individual case studies or small
groups of populations by, for example, Sarro et al. (2014); Bouy et al.
(2015); Herczeg & Hillenbrand (2015); Olivares et al. (2018b); Miret-Roig et
al. (2019); Olivares et al. (2019); Li et al. (2021); Olivares et al. (2021,
2023). We now want to take this effort a step further and create an archive of
empirical isochrones for a large sample of open clusters in the solar
neighborhood, using purely observational data combined with a homogeneous
extraction method based on statistical analysis and machine learning tools. By
sacrificing the ability of absolute age quantification, we can decouple
isochrones from their strong associated model dependencies and still have a
relative age determination method that is entirely based on observational
data. This approach is also independent of any assumptions regarding
extinction or metallicity. By comparing new and known clusters, our archive
can provide insights into cluster evolution and the properties of large-scale
structures in the galaxy. Extracting and analyzing empirical isochrones may
lead to the discovery of new features in their shapes, particularly regarding
the lower main sequences, and provide calibration data for stellar evolution
models, improving their quality and broadening our understanding of stellar
evolution.

##  3 Data

In the following paragraphs, we outline the surveys and cluster member lists
that we consulted for creating the empirical isochrone archive. As their
specific membership determination method often influences our choice of
applied quality cuts, we briefly describe the applied techniques used by the
respective data source.

###  3.1 Survey data

The observational information for the archive clusters is collected almost
exclusively from the data products of the _Gaia_ mission (Gaia Collaboration
et al. 2016). Specifically, the extraction algorithm was built on input from
the early third data release (EDR3, Gaia Collaboration et al. 2021), as well
as the full DR3 (Gaia Collaboration et al. 2022), which contain identical
astrometric and photometric measurements. To facilitate comparisons to pre-DR3
works, we additionally include isochrones based on _Gaia_ DR2 (Gaia
Collaboration et al. 2018b) in our archive.444All CMDs using _Gaia_ passbands
shown throughout this work depict _Gaia_ DR3 data unless explicitly indicated
otherwise.

We chose the _Gaia_ survey as our primary data source mainly due to its
unprecedented size and precision. Its latest data releases include five-
parameter astrometric solutions for 1.2 and 1.5 billion sources, respectively.
Profiting from a homogeneous source of photometric and astrometric
measurements for all archive clusters, we greatly reduce the need for
gathering data across multiple sources (Cantat-Gaudin et al. 2018), which has
been reported to cause discrepancies in the results of cluster analyzes of
different authors (e.g., Netopil et al. 2015). Furthermore, the low parallax
uncertainties of $0.1$ $\mathrm{m}\mathrm{a}\mathrm{s}$ for sources with G
$\lesssim$ 18 mag, and the sub-millimagnitude photometric uncertainties in the
color measurements ( $\leq$ 13 mag) provide the best possible conditions for
creating sharp color-magnitude diagrams of star clusters. Finally, due to its
homogeneous nature and vast source catalog, the _Gaia_ releases have been
established among the most commonly used optical surveys for performing
cluster membership and isochronal age determinations.

Alongside the _Gaia_ data, we also draw on data from the Two Micron All Sky
Survey (Skrutskie et al. 2006, 2MASS) $J,H,K_{S}$ photometric system, as well
as the Panoramic Survey Telescope And Rapid Response System (Chambers et al.
2016, Pan-STARRs) $u,g,r,i,z,y$ bands for the two nearby clusters IC 4665 and
the Pleiades. We do this to test the code on photometric systems other than
the one it was designed on and to study empirical lower main sequences down to
the mass level of brown dwarfs.

The required parameters for the extraction of empirical isochrones are
constrained to apparent magnitude or flux measurements in at least two
passbands for calculating the color indices, and stellar distances, to infer
the absolute magnitudes via the distance modulus. The distances are computed
via an inversion of the parallax, rather than for example using a Bayesian
method, such as the one proposed by Bailer-Jones (2015), or the systematic
parallax offset described in Lindegren et al. (2018). To justify this choice,
we analyzed the fractional parallax errors of all datasets where we performed
the distance estimation ourselves, meaning the _Gaia_ DR2 and DR3 source
catalogs. We found that only for 0.03 % of all DR3 sources and 0.42 % of all
DR2 sources do those errors exceed 10 %. Therefore, distance calculation via
parallax inversion is a viable approach. However, we note that in case of
using datasets with typical parallax errors larger than 10 %, we recommend an
alternative way of distance estimation when calculating absolute magnitudes.

###  3.2 Cluster data and membership lists

Table 1: Overview of the catalog contents, mean cluster parameters, and the
applied quality cuts.

Member | Cluster(s) | Mean distance | Member stars | Probability | Age ref. | Further cleaning
---|---|---|---|---|---|---
selection |  | (pc) | $N_{*}$ | or score |  |
_Catalogs_
(1) | 64 | $135-498$ | $\leavevmode\nobreak\ 96-1749$ | $p\geq 0.5$ | (a), (b) | Q1
(2) | 10 | $137-402$ | $321-1828$ | no | (a) | Q1
(3) | 16 | 1$06-177$ | $150-1196$ |  stability $>6$ | (c) | Q2; $G_{\mathrm{BP,err}}<0.05$ mag
_Single clusters_
(4) | Hyades | 67 | 972 | no | (a) | Q1
(5) | Coma Berenices | 86 | 181 | no | (d) | Q1
(6) | Meingast 1 | 135 | 1366 |  stability $\geq 0.24$ | (a) | Q1
(7) | IC 4665 (DANCe) | 356 | 819 |  $p_{\text{\emph{Gaia}}}>0.7$ or $p_{\text{DANCe}}>0.5$ | (a) |  $i>13$ mag; $r>9$ mag
(8) | Pleiades (DANCe) | 136 | 1483 | $p_{t}\geq 0.84$ | (a) |  $e_{\mathrm{i}}<0.3$ mag, $z>0$ mag

_Member selections._ (1) Cantat-Gaudin & Anders (2020); (2) Meingast et al.
(2021); (3) Ratzenböck et al. (2022); (4) Meingast, priv. comm.; (5) Fürnkranz
et al. (2019); (6) Ratzenböck et al. (2020); (7) Miret-Roig et al. (2019); (8)
Olivares et al. (2018a); _Ages._ (a) Cantat-Gaudin et al. (2020); (b) Dias et
al. (2021); (c) Ratzenböck et al. (2023); (d) Curtis et al. (2019) 555All
distances rounded to integer values. Numbers and quality criteria referring to
DR3 data.
Quality filter 1 (Q1): $\varpi>0$, RUWE${}_{\mathrm{DR3}}\leq 1.4$.
Quality filter 2 (Q2): $\varpi>0$, RUWE${}_{\mathrm{DR3}}\leq 1.4$, fidelty_v2
$>0.5$, $G_{\mathrm{err}}<0.007$ mag, $G_{\mathrm{RP,err}}<0.03$ mag
(Ratzenböck et al. 2022, 2023).

We collect cluster memberships from various sources found in contemporary
literature, under the consideration of two global quality cuts: First, the
mean cluster distance is limited to $d\leq$
$500\text{\,}\mathrm{p}\mathrm{c}$, and secondly, the number of high fidelity
cluster members is required to be $N_{*}\geq 100$. By “high fidelity” we group
the probability or stability scores calculated by different methods, as can be
seen in the fifth column of Table 1.

The first cut is justified from an observational point of view, as sources
grow increasingly faint as a function of the cluster distance, resulting in an
ever-increasing under-representation of low-mass sources in the CMDs of more
distant clusters. At the same time, _Gaia_ parallax errors increase
disproportionately when moving out to farther distances. For instance, Cantat-
Gaudin & Anders (2020) state that the uncertainty in the proper motion
measurements of DR2 starts to contribute significantly to the scattering of
the cluster members already at distances around $\sim$
$500\text{\,}\mathrm{p}\mathrm{c}$; hence we pick this distance as our cutoff
value.

The member number limit originates from the extraction method’s perspective:
Empirical isochrones should ideally represent a significant portion of a
cluster’s stellar initial mass function (IMF) and trace its population over a
wide dynamical range. As clusters age, low-mass members may be dynamically
ejected, leading to a depletion in the low-mass region of the main sequence.
However, while this mass loss and its treatment in different membership
determination methods are noteworthy, the general abundance of low-mass stars,
owing to the shape of the IMF, ensures a well-populated lower main sequence.
Consequently, the proposed empirical isochrone extraction is not significantly
affected by a slight change of the mass distribution in the low-mass regime.

Another motivation of the member number cut stems form the fact that the
extractions are performed with statistical methods, which require a
sufficiently large set of data points for each cluster to yield reliable
results. From heuristic trial-and-error experiments, we determined that a
member list of at least 100 sources provides adequate CMD populations in the
overwhelming majority of the cases.666After the cross-matching and data
cleaning, a cluster occasionally ends up with less than 100 members. These
variations of a maximum $\pm 5\leavevmode\nobreak\ \%$ are allowed within our
cut. Depending on the sharpness of the CMD of an observed population, one
could also relax this criterion, but it would come with the stipulation of
having to inspect each case individually manually.

We accumulate 100 membership lists corresponding to 88 individual star
clusters. The discrepancy between the number of memberships and actual
clusters in the archive arises from a twofold overlap: On the one hand, the
selection of 10 clusters from Meingast et al. (2021) corresponds to a subset
of the sample collected from Cantat-Gaudin & Anders (2020). On the other hand,
for the two clusters IC 4665 and the Pleiades, we additionally collect data
from the Dynamical Analysis of Nearby Clusters (DANCe) survey (Olivares et al.
2018a; Miret-Roig et al. 2019) along with the _Gaia_ observations. Concerning
the member selections based on the _Gaia_ passbands, only the selection of
Ratzenböck et al. (2022) was performed on DR3 data, whereas all others used
DR2 data. Since the third data release generally comprises more accurate
measurements, we cross-match the DR2 selections with their DR3 counterparts.

Due to quality concerns regarding the empirical isochrones that the algorithm
calculated in isolated cases, we excluded five clusters from our original
sample (details in Sect. 3.3), bringing the number of archive clusters to 83
in total. A summary of the cuts and mean cluster parameters of our final
selection, divided into their respective source catalogs, can be found in
Table 1. The following provides a brief overview of the literature sources and
a short description of their membership determination processes.

####  3.2.1 Cantat-Gaudin & Anders (2020) selection

Our primary source is the extensive open cluster catalog by Cantat-Gaudin &
Anders (2020), who created membership lists for 1481 clusters in the Milky
Way. After applying the global quality cuts, with high fidelity sources being
defined as $p\geq 0.5$ on the probability score determined by the authors, we
obtain membership lists for 67 open clusters, hereafter referred to as Catalog
I. Of these, 59 were first determined in Cantat-Gaudin et al. (2018) with the
unsupervised, two-step UPMASK procedure (Unsupervised Photometric Membership
Assignment in Stellar Clusters, Krone-Martins & Moitinho 2014). Their method
consisted of clustering stars based on astrometric proper motion and parallax
measurements ($\varpi$, $\mu_{\alpha_{*}}$, $\mu_{\delta}$) of _Gaia_ DR2, and
subsequently validating each cluster by comparing its density to that of a
random distribution. In doing so repeatedly with values drawn from the
astrometric parameter distributions, a membership probability was determined
for each star in their input sample (Cantat-Gaudin et al. 2018). Catalog I
further includes 5 UPK clusters (Sim et al. 2019), for which Cantat-Gaudin &
Anders (2020) also provide membership lists calculated with UPMASK. Finally, 3
of the 67 groups are UBC clusters found by Castro-Ginard et al. (2018), who
used DBSCAN (Ester et al. 1996) in combination with artificial neural networks
on the 5D astrometric parameter space of _Gaia_ DR2 data to perform their
selection. It should be noted that Castro-Ginard et al. (2018) did not compute
membership probabilities for their clusters, and thus the probability column
provided in Cantat-Gaudin & Anders (2020) equates to one for all UBC clusters.

####  3.2.2 Meingast et al. (2021) selection

For a subset of Catalog I, namely the ten open clusters $\alpha$ Per (Melotte
20), Blanco 1, IC 2602, IC 2391, NGC 2451A, NGC 2516, NGC 2547, NGC 7092,
Platais 9, and the Pleiades we also collect alternative membership lists from
Meingast et al. (2021), hereafter referred to as Catalog II. These authors
determined cluster membership based on deprojected proper motion measurements,
which enabled them to not only retrieve members in the cluster cores but also
in sub-field-density stellar coronae that can extend several 100 pc from the
central region. As a result, comparatively more members were detected,
especially down to the low-mass end of the main sequence. With their
clustering method, no membership probability was assigned to the members of
the clusters.

We use this subset of open clusters as main sample for testing and refining
the extraction algorithm. Furthermore, as the lower main sequences of their
clusters are generally better populated than those of the Catalog I clusters,
they are well suited for studying the information gain via empirical
isochrones in the very faint regions (see Sect. 6.2.2).

####  3.2.3 Ratzenböck et al. (2022) selection

To extend the empirical isochrone archive toward young clusters, we consult
the work of Ratzenböck et al. (2022), who recently extracted 37 co-moving star
clusters in the nearby Sco-Cen star-forming region. To determine the cluster
memberships, Ratzenböck et al. (2022) employed the Significance Mode Analysis
(SiGMA) tool to _Gaia_ DR3 data in a box around the molecular cloud complex.
SiGMA works by identifying modal regions in the 5D astrometric parameter space
of the observational data. It first estimates the data density and performs a
graph-based gradient hill climb to identify local density peaks. The multiple
resulting preliminary clusters are then iteratively merged using a modality
test to measure the significance of the separating density dip. Instead of
membership probabilities, the algorithm assigns each cluster member a
stability value between zero and 100, though a given cluster may have a
maximum stability of less than 100. We choose a generous threshold of
stability $>6$ for high fidelity members. Furthermore, we apply the quality
cuts listed in Ratzenböck et al. (2023) (see Table 1) and set a stricter
quality cut in the blue passband error $G_{\mathrm{BP,err}}<0.05$ mag to
remove outliers from the lower left region of the CMD. As a result, we retain
18 of their determined clusters, which we hereafter refer to as Catalog III.

Figure 3: Spatial distribution and relative sizes of the archive clusters,
determined by the different source catalogs. The cluster positions included in
Catalog II are not displayed to avoid duplication of Catalog I sources. The
marker colors correspond to the estimated cluster ages from literature (Table
1).

####  3.2.4 Hyades and Coma Ber clusters

To complement the clusters of Catalog I, we collect membership lists for the
Hyades (Melotte 25) and the Coma Berenices (Melotte 111) cluster. The two
well-known OCs are absent from the works of Cantat-Gaudin et al. (2018);
Cantat-Gaudin & Anders (2020), as their proximity makes those clusters
incompatible with their membership selection method. Regarding Coma Ber, we
use the membership selection published by Fürnkranz et al. (2019), who
performed a wavelet decomposition to single out stellar overdensities in
velocity space and subsequently employed DBSCAN to identify cluster members
without the use of a probability or score measure. For the Hyades, the
clustering method is the same that was employed for the clusters of Catalog II
(Meingast, priv. comm.).

####  3.2.5 Meingast 1 stellar stream

The Meingast 1 stream (Meingast et al. 2019; Curtis et al. 2019; Ratzenböck et
al. 2020; Röser & Schilbach 2020, also referred to as Pisces-Eridanus stream)
is not a classical open cluster, but rather a large, spread-out stellar group
of co-moving objects, which traverse the solar neighborhood at roughly a
distance of $100\text{\,}\mathrm{p}\mathrm{c}$. At a length of more than
$400\text{\,}\mathrm{p}\mathrm{c}$ and width of
$50\text{\,}\mathrm{p}\mathrm{c}$, it is speculated to be the remnant of a
tidally disrupted cluster or OB association. The stream was initially
identified via a wavelet decomposition of the 3D velocity space of stars in
the solar neighborhood ($d<$ $300\text{\,}\mathrm{p}\mathrm{c}$) by Meingast
et al. (2019). In a later work, Ratzenböck et al. (2020) revisited the object
and determined new memberships based on the five-dimensional DR2 position and
velocity space ($\alpha$, $\delta$, $\varpi$, $\mu_{\alpha}$, $\mu_{\delta}$)
clustering of a supervised one-class support vector machine (OCSVM, Cortes &
Vapnik 1995). Similar to the clusters of Catalog III, possible member stars
were assigned a stability parameter, albeit normalized and provided as
percentages for the stream. Of the two filter criteria proposed by Ratzenböck
et al. (2020) (stability $>4$ % or stability $\geq 24$ %), we apply the
stricter one.

As the stream members are still thought to be of a common origin, displaying a
narrow main sequence very similar to the Pleiades (Ratzenböck et al. 2020), we
also choose to test our isochrone extraction algorithm on this object. We
further compare the empirical isochrones of the stream with those of OCs of a
similar evolutionary state in Sect. 6.2.2 to analyze the impact of empirical
isochrones on the lower main sequence in the isochrone blindspot region.

####  3.2.6 DANCe data for IC 4665 and the Pleiades

The data for the clusters IC 4665 and the Pleiades was collected as a part of
the DANCe survey (Bouy et al. 2013; Olivares et al. 2018a; Miret-Roig et al.
2019). As such, the measurements for the clusters were gathered from different
observatories and instruments and contain observations in near-infrared
passbands. Since their data do not include parallaxes, these values need to be
substituted. For IC 4665, Miret-Roig et al. (2019) already calculated a median
parallax value from the portion of cluster members visible in _Gaia_ DR2.
Concerning the Pleiades, we adapt the distance from the cluster parallax
denoted in Cantat-Gaudin & Anders (2020). Regarding the membership
determination methods, Miret-Roig et al. (2019) used an algorithm that models
a field and a cluster component, respectively, based on an initial member list
(Sarro et al. 2014; Olivares et al. 2019). Their approach is of Bayesian
nature, and a membership probability is provided for the clustered sources.
The membership determination for the Pleiades was undertaken with a similar,
but more computationally expensive, Bayesian-based technique, consisting of a
generative mixture model and a Markov Chain Monte Carlo method for calculating
the posterior distributions of the parameters. A probability score for cluster
members was created by adding the Bayesian membership probability and its
sensitivity to the cluster parameter for each star in their sample. Their
probability threshold was defined at $p_{t}=0.84$ (see Table 1), which we
adopt for the high fidelity sources.

Including the two clusters in our selection offers a convenient way to
accomplish studying the shapes and added value of empirical isochrones
regarding the faint, low-mass end of the main sequence of star clusters.
Additionally, it permits exploring the extension of the method toward
different photometric systems.

###  3.3 Final cluster selection and general archive properties

After applying our method to all 88 clusters described in the previous section
and visually inspecting the results and the respective cluster CMDs, we had to
perform a final cut on our archive selection. This is because, for a few
select sources, no physically meaningful isochrones could be extracted from
their CMDs due to either a severe distribution scattering or concerns
regarding their member selections. The concerned objects include the three
clusters L1641 S, IC 348, and RSG 7 from Catalog I, as well as the two groups
$\rho$ Oph/L1688, and Lupus $1-4$ from Catalog III.

Three of the excluded clusters are estimated to be very young (L1641 S $\sim$
$3\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$; Da Rio et al. 2017, $\rho$
Oph/L1688 $\sim$ $3.1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$, and Lupus $1-4$
$\sim$ $4.5\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$; Ratzenböck et al. 2023),
which is likely the reason for their scattered CMDs.

On the other hand, with an estimated age of $\sim$
$12\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ (Cantat-Gaudin et al. 2020), the
cluster IC 348 does not fall below the apparent lower age threshold of the
archive. However, its estimated extinction of $A_{V}=1.91-2.3$ mag exceeds
that of the second highest reported extinction (UBC 17a, $A_{V}=0.8$ mag) by
almost a factor of two (references from Table 1). The region was for example
studied by Olivares et al. (2023), who further determined that the extinction
varies significantly over the dynamical range of the cluster CMD and discussed
its extinction law. We note that along with $A_{V}$ and $R_{V}$ distributions,
the authors also provide empirical, extinction-free isochrones derived with a
similar method for IC 348 and five other Perseus groups in their work.

In contrast to the other excluded clusters, RSG 7 is not characterized by any
extreme values concerning either age or extinction. Still, following our
quality cuts, its CMD displays two distinct populations along the main
sequence, directly contrasting the sharp, well-defined age sequence one would
expect for a young, bound cluster. As both subbranches are comparable in
member numbers, it is impossible to distinguish between contaminants and
actual cluster members without further investigation. Therefore, we exclude
this cluster.

Figure 3 displays the spatial distribution of the final 83 clusters selected
for the archive in a heliocentric Cartesian coordinate system. The number of
cluster members, determined by the respective source catalog, is reflected in
the relative size of the position markers. The clusters are color-coded
according to the estimated ages found in the literature. The references for
the cluster ages are provided in the rightmost column of Table 1. For the
cases of Catalog I, we preferably use the ages from Cantat-Gaudin et al.
(2020) and only consult the work of Dias et al. (2021) as a secondary source
for clusters not covered by the former. A histogram showing the age
distribution of the archive clusters is displayed in Fig. 4. As can be seen,
the bulk of the archive clusters are younger than
$100\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$, and only very few clusters are
older than one Gyr, whereas the lower age limit is around 7 Myr. The displayed
age caps indicate a reliable age range for a representative cluster catalog
and correspond to natural limits in open cluster observations. After all, very
young clusters ($\lesssim$ $5\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$) are at
least partly obscured and reddened by their parental dust cloud, rendering
them invisible to _Gaia_ or badly scattered, as discussed above. Very old
clusters ($\gtrsim 1$ Gyr), on the other hand, are naturally rare due to
infant mortality (Lada & Lada 2003, and references therein).

Figure 4: Histogram depicting the logarithmic age distribution of the archive
clusters, according to the literature values listed in Table 1. Due to
visualization purposes, the ordinate is scaled logarithmically. The ages are
binned using Knuth’s rule, as is implemented in astropy (Astropy Collaboration
et al. 2013, 2018).

##  4 Methods

In search of an isochrone extraction routine, we tested several methods, some
of which have their roots in data analysis, while others originate from the
research field of image processing. In addition, we investigated various
preprocessing strategies. After evaluating their respective performance on the
clusters of Catalog II, we decided on a combination of Support Vector
Regression (SVR) and Principal Component Analysis (PCA) as the final
procedure. Not only did this combination yield the best results, but it also
conforms to an important applicability aspect, as it is very adaptive toward
the many different CMD shapes in our sample. Thus, it can be easily applied to
CMD distributions reflecting various evolutionary states of the corresponding
clusters and diagrams created from different passband combinations or
photometric systems, without requiring extensive manual parameter tuning. We
employ PCA as preprocessing step to provide continued functional dependency of
the cluster data, particularly concerning evolved clusters. SVR is then used
to perform the actual extraction of the empirical isochrones. While other
authors have also started providing empirical isochrones for individual case
studies of open clusters, to the best of our knowledge, we are the first to
design an approach utilizing SVR and to derive empirical curves for a large,
age-representative sample of clusters in a homogeneous manner. For comparison,
Bouy et al. (2015); Miret-Roig et al. (2019); Olivares et al. (2019) for
instance employ a version of the Principal Curves technique first described by
Hastie & Stuetzle (1989), which we also tested but ultimately discarded due to
the loss of some key characteristics of the isochrones. Both PCA and SVR are
implemented in the Python library scikit-learn (Pedregosa et al. 2011) in
their nonlinear extensions, which we use for our algorithm.

###  4.1 Requirements and possible pitfalls

Despite the major benefit of model independence, one should be aware of some
potential pitfalls when relying on observational data to produce isochrones.
For instance, a low signal-to-noise ratio in the measurements can lead to a
large scatter in the star positions in the CMD. This makes it difficult to
derive the finer characteristics of the isochronal curve empirically. Another
common occurrence in observational data concerns outliers. In cluster CMDs,
they often appear in the form of unresolved main sequence binaries (see Sect.
2). Unresolved binary pairs of main sequence stars with evolved objects such
as red giants or white dwarfs can induce further scatter both above and to the
right, as well as below and to the left of the main sequence in a CMD (Gaia
Collaboration et al. 2018a). Field stars can cause additional contamination of
the member selection, either assuming random positions in the CMD or forming a
distinct sequence (e.g., RSG 8 in Fig. 15). Furthermore, all observational
data are associated with measurement errors, which in the case of _Gaia_ DR3
are readily available and should be considered in the isochrone extraction. A
final caveat arises from the fact that the extraction technique is based on
data analysis, which generally requires a functional dependency between the
variables one is interested in. Common isochrone characteristics such as the
main sequence turn-off, white dwarf regions, or even a relatively vertical
distribution of upper main sequence stars directly counteract such a
functional dependency in the original CMD parameter space and must be
accounted for in the algorithm design.

###  4.2 The algorithm

To counteract the described possible caveats, the cluster data need to be
preprocessed in a way that preserves information while simultaneously ensuring
a continued functional relationship between the input variables. Furthermore,
the isochrone extraction needs to be robust against outliers assuming various
forms and needs to be able to incorporate measurement errors. Hence, the
workflow of creating a robust empirical isochrone from observational data with
our developed method consists of the following main steps:

  1. 1.

Preprocessing via PCA transformation of the cluster data.

  2. 2.

SVR hyperparameter tuning via gridsearch and cross-validation.

  3. 3.

Extraction of a preliminary empirical isochrone with the tuned SVR model and
the weights calculated from measurement errors.

  4. 4.

Generation of a set of resampled isochrones via iterative bootstrapping
($n_{\text{boot}}=1000$).

  5. 5.

Calculation of the median, the 5${}^{\text{th}}$ and 95${}^{\text{th}}$
percentile values from the collection of resampled and reverse-transformed
isochrones.

In the following paragraphs, each step of the workflow is discussed in more
detail.

####  4.2.1 Preprocessing: Principal component analysis

Principal component analysis is an adaptive data analysis technique (Jolliffe
& Cadima 2016), whose variables are always defined by the specific input
dataset instead of a priori. Thus, it is well suited for an isochrone
extraction algorithm that needs to perform on a versatile input cluster CMDs.
Briefly put, PCA transforms the input parameters, in this case the color index
and absolute magnitude of a cluster CMD, into new, uncorrelated variables that
successively maximize the variance in the data. It can be shown that this
maximization can be reached by defining the new variables as a set of mutually
orthogonal linear combinations of the input variables, together with the
eigenvectors corresponding to the largest eigenvalues of the data. These
combinations are commonly referred to as “principal components.” For an
introduction to the basic principle of the method, we refer readers to works
such as Tipping & Bishop (1999); Jolliffe & Cadima (2016). A detailed
description of the application of the technique to the problem at hand can be
found in Appendix C.

Figure 5: Illustration of a PCA transformation of CMD data on the example of
the Pleiades cluster. _Left:_ Extent and orientation of the principal
components. The first component (_green_) is much larger than the second one
(_magenta_), reflecting the different orders of magnitude of variance in the
observation matrix. The directions of the PCA space axes are denoted by black
arrows, with the second one magnified by a factor of 50 for better visibility.
_Right:_ Transformed cluster data in PCA space, with the new axes pointing in
the direction of the two principal components.

We use PCA to rotate the reference frame into the direction of the principal
components spanned by the new variables (see Eq. 1). The purpose of this data
transformation is illustrated in Fig. 5 on the example of the Pleiades
cluster. The left panel shows a regular CMD, overlaid by two colored bars
indicating the extent of the principal components, which have been scaled by
their corresponding variance. The arrows represent the direction of the
eigenvector associated with the respective component and therefore span the
axes of the PCA space. Due to the source distribution in most CMD
configurations, which mostly comprises stars populating the diagonally or
vertically inclined main sequence region, the variance in the direction of the
absolute magnitude axis is much larger than the one for the color axis. This
physical fact is accurately represented in the ratio of the principal
components; it can be seen that the second principal component is much smaller
than the first one, and the superimposed arrow denoting the direction of its
eigenvector has been magnified by a factor of 50 for better visibility. The
PCA-transformed cluster data can be seen in the right panel of Fig. 5. A look
at the $y-$axis of the PCA space shows that the absolute data scatter is much
smaller in that direction than for the CMD. In preprocessing the cluster data
with PCA, we circumvent two of the pitfalls of observational data discussed in
Sect. 4.1: Firstly, projecting the data into PCA space concentrates the
variance on the first component instead of distributing it between the two
input variables, which is beneficial for the isochrone extraction. Secondly,
the functional dependency between the two variables is preserved for all
shapes and evolutionary stages of cluster CMDs.

####  4.2.2 Curve extraction: Hyperparameter tuning and weighted support
vector regression

The central component of the algorithm, meaning the method governing the
generation of empirical isochrones, is based on the machine learning technique
of support vector regression. The sparse, memory-efficient approach originates
in statistical learning theory (Smola & Schölkopf 2004) and was, for example
described by Bishop (2006), to whom interested readers are referred for a more
detailed description. The procedure’s working principles are outlined in
Appendix D, based on its implementation in sklearn.

The general benefits of the method include its straightforward implementation,
low computational cost, and small number of hyperparameters, meaning variables
not subjected to analytical optimization. Instead, they can be defined
somewhat arbitrarily and impact the resulting regression to varying degrees.
The hyperparameters of SVR originate at different points of the formal
derivation of the method (see Appendix D for a detailed discussion). The first
two, epsilon and C, are defined via the characteristic loss function of the
regression (Eqs. 8$-$11). The former determines the number of support vectors,
meaning the size of the subset of training data at which the kernel function
is evaluated. By adjusting its value we can enhance the method’s sensitivity
for outliers, particularly regarding the unresolved binary sequence. The
latter parameter, also called the penalty, determines the leniency toward
considering points far away from the main distribution in the regression.
Lastly, for our algorithm we use a radial basis function as kernel, adding a
final hyperparameter called gamma, which is connected to the variance of the
input data. Due to the many possible CMD shapes, the hyperparameters must be
tuned once for a given CMD variation of a cluster at the beginning of a
regression. However, they do not need to be re-established after the isochrone
has already been calculated for a specific CMD or group of similar CMDs. We
tune the hyperparameters on a predefined parameter grid using a gridsearch and
5-fold cross-validation within the sklearn package. For a detailed analysis of
the influence of each hyperparameter and a description of the tuning process,
readers are referred to Appendix D.2.

Once the hyperparameters for a given CMD configuration have been determined,
the tuned SVR model is fitted to the entire transformed cluster data,
corresponding to the whole set of observations. This is done as we are not
predicting further points but instead want to generate a regression curve
representing the observed data as accurately as possible. Using the trained
model to predict the response variable, meaning the PCA Y data, yields a
regression curve that can then be easily transformed back into the original
parameter space, where it corresponds to an empirical isochrone, akin to the
data transformation shown in Fig. 5. A benefit of the SVR implementation of
scikit-learn (Pedregosa et al. 2011) is its ability to incorporate a weight
parameter for each data point, which is subsequently used to multiply the
penalty parameter C. Large weights urge the SVR model to include a data point
in the regression, whereas small weights allow more lenience. Hence, we can
solve the problem of incorporating measurement errors, if they are available,
by turning them into one-dimensional weights for each data point. To this end,
we convert the photometric and astrometric measurement uncertainties into
scalar values using the root sum of squares (RSS), as well as Gaussian error
propagation (details in Appendix D.3).

####  4.2.3 Empirical isochrones and uncertainty regions

Figure 6: Illustration of the bootstrapping and resampling process. The
procedure is shown for the well-behaved cluster NGC 2422 (_top row_), and the
more scattered beta Sco group (_bottom row_). The _orange lines_ in the left
column of the figure indicate 100 resampled SVR isochrones created via
bootstrapping of the PCA cluster data. The center column depicts the same
lines projected into the respective cluster CMD. The right column shows the
median (_red line_), 5${}^{\text{th}}$ and 95${}^{\text{th}}$ percentile
boundaries (_black lines_) calculated from the resampled set.

Figure 7: Summary plot of all _Gaia_ empirical isochrones calculated for the
83 archive clusters, shown for the three different color indices
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$,
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}$, and
$\mathrm{G}-\mathrm{G}_{\mathrm{RP}}$. The color coding of the isochrones
corresponds to the estimated cluster age found in the literature (see Table 1
for details). An interactive version of the figure is also available here.

The last two steps of the algorithm workflow are designed to create a more
reliable empirical isochrone than can be achieved by a single regression.
Therefore, we bootstrap the PCA-transformed cluster data to create a slightly
different input array for training a new SVR model with the predetermined
hyperparameters. To get a robust result, we create $n_{\text{boot}}=1000$
bootstrapped datasets for training SVR models and calculate the regression
curves in PCA space by predicting the original PCA dataset. We then use the
set of resampled curves to determine the empirical isochrone as the median of
the values for each color index of the original data. We further compute the
5${}^{\text{th}}$ and 95${}^{\text{th}}$ percentile values, respectively. They
correspond to an uncertainty region around the empirical isochrone and
represent its reliability or trustworthiness locally along the data
distribution in the CMD. The bootstrapping and the median or percentile array
calculations are undertaken in PCA space before the results are inversely
transformed into the original CMD space using the transformation matrix
defined for the original cluster data.

The process is visualized in detail for two example clusters in Fig. 6, using
only $n_{\mathrm{boot}}=100$ bootstrapped curves for better visibility. For
well-sampled clusters, such as the one in the _top row_ of the figure, the
individual curves are generally in excellent agreement and largely
overlapping. As a result, the uncertainty regions are very narrow over the
whole populated sequence. On the other hand, if the cluster distribution is
more scattered or exhibits sparsely populated regions, as in the case of the
young beta Sco group, those areas are characterized by larger uncertainty
bounds, which tend to broaden toward the upper and lower limits of the
empirical isochrone.

##  5 Results

We derived empirical isochrones for the three CMD color combinations of _Gaia_
DR2 and DR3 for all 83 clusters in the final source catalog (Sect. 3.3). The
tables containing the results and the calculated uncertainty bounds are
available online; an example of their contents can be seen in Appendix B.

The results based on the third _Gaia_ data release are displayed in Fig. 7 and
color-coded according to the estimated cluster age found in the literature
(Table 1), akin to the spatial distribution shown in Fig. 3. An interactive
version of the summary figure, in which isochrones can be selected arbitrarily
and studied in more detail, or their evolution may be visualized by
successively adding them to the diagram using a slider corresponding to their
estimated ages, is also available online and can be accessed via the link in
the figure caption. Besides the summary plot, figures 15 and 16 provide a
separate view of the results for each cluster, using a
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ color-magnitude diagram
and the corresponding empirical isochrone and uncertainty bounds. The panels
displaying the clusters are sorted in ascending order of their estimated
literature ages.

From the individual cluster CMDs, it is apparent that the SVR-generated
isochrones generally trace the respective stellar distributions very well.
They even preserve specific details, such as the knee-like bend in the upper
main sequence ($\mathrm{M}_{\mathrm{G}}\sim 4$ mag,
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}\sim 1-1.5$ mag)
characteristic for young stellar populations, in many cases, even though this
part of the CMD is often sparsely populated. Regarding the lower main sequence
region, the empirical isochrones pass through the most densely populated part
of the stellar distribution, which we assume to be the most reliable indicator
of the true distribution of the cluster stars, that they would assume without
the scatter induced by observational errors. Another characteristic captured
by the empirical isochrones is the arc toward the upper right of the CMD of
evolved clusters. This feature is reminiscent of a main-sequence turn-off,
albeit not as intricate as ones available from model isochrones due to the
lack of observable sources at this evolutionary stage. Likewise, the empirical
indicator of a main sequence turn-off is only present when our cluster
selection includes observed members in the red giant region. This is why some
evolved clusters with similar literature ages do not exhibit a turn-off in the
empirical isochrones while those in adjacent panels of Fig. 16 do. This
circumstance should also be considered when looking at the number of main
sequence turn-offs visible in the summary figure.

Overall, Fig. 7 displays a generally good agreement between the distribution
of the empirical isochrones and the cluster ages found in literature, as
indicated by the mostly smooth color gradient visible between the isochrones
depicting clusters that were determined to be young and those portraying ones
that have been estimated to be older. Looking at the evolved clusters for
which we do recover a main sequence turn-off indicator, we also report a good
agreement between the model-determined ages and the empirically determined
ones based on the evolution of the turn-off point positions, although a few
disagreements are also apparent. For instance, the Stock 2 cluster has a lower
empirical isochrone turn-off, which indicates a more evolved state than its
literature age suggests. On the opposite end of the age range, it can be
observed that the empirical isochrones corresponding to the youngest clusters
in our sample successfully and unanimously capture the relative shift of the
lower main sequence toward the upper right of the CMD.

Overall, the good agreement in the comparison between the literature cluster
ages determined via model fitting and the purely observation-based empirical
isochrones provides an important mutual verification for both approaches.
Apart from possible systematic errors, model-dependent cluster ages from
isochrone fitting do indeed often represent the stellar evolution of clusters
very well. On the other hand, our archive of empirical isochrones proves to
cover an extensive range of evolutionary phases of clusters. Both facts are
apparent from the correlation between the cluster age from literature,
indicated by their color, and the locations of the corresponding empirical
isochrones in Fig. 7. Regardless, this strong correlation breaks down for
individual clusters.

###  5.1 Validation using Sco-Cen clusters

To assess the practical feasibility of using the empirical isochrone archive
as an age-scaling ladder, we validate our results against the clusters of
Catalog III. The dataset offers the best testing conditions of our source
catalogs, as its populations all originate from the same star-forming region,
and their ages have been determined homogeneously. Additionally, of the 16
considered groups, only beta Sco has a reported mean extinction greater than
zero ($A_{V}=0.5$ mag). Thus, the validation is neither biased against
different age estimation methods nor against differing general parameters
between the various groups. Moreover, with estimated ages between circa 7 and
20 Myr, the CMD distributions of the cluster isochrones are expected to show a
significant shift of the lower main sequence toward the upper right, which
should make them clearly distinguishable from one another. At the same time,
their low absolute age spread means that they would provide a very fine grid
for a possible age-scaling of young, unknown stellar populations.

Figure 8: Summary plot of empirical isochrones derived for the 16 clusters out
of the 34 groups discovered in the Sco-Cen complex (Ratzenböck et al. 2022),
which are in the empirical isochrone archive. The lines are color-coded in
accordance with the estimated ages published in Ratzenböck et al. (2023). The
age sequence constitutes an almost perfect fit to the empirical isochrone
shapes and positions within the CMD.

Figure 8 displays the calculated empirical isochrones for the Sco-Cen groups,
with the different colors again indicating their respective estimated ages
(Ratzenböck et al. 2023). As for Fig. 7, a strong correlation between the
locations, meaning the evolutionary phases, of the empirical isochrones and
the associated ages of the populations is visible for the entire sample. By
building up the plot sequentially, we found that in the region between $\sim
4$ and $11$ G magnitudes, where the uncertainty bounds of the empirical
isochrones are the narrowest, there is even a perfect agreement between the
isochrone position and age sequence. The trend also holds true in the upper
main sequence ($<2$ mag) for most of the groups, except for the sigma Sco and
the Antares clusters. A visual inspection of the CMDs of the two populations
showed that their distributions exhibit a higher apparent scatter in their
upper main sequences than the other Sco-Cen groups. As a result, the algorithm
cannot derive a physically meaningful upper main sequence region for the two
clusters.

From the favorable correlation, we conclude that the Sco-Cen empirical
isochrones may indeed be used as an age-scaling grid or ladder for determining
relative ages of other, unknown populations with similar ages. Given the
excellent verification between literature and empirical isochrone positions in
Fig. 7, the principle can be extended to the whole archive, though extra care
needs to be taken when comparing between clusters with ages determined by
different authors, as those estimates might be biased again by their
respective isochrone fitting or model choices.

###  5.2 Isochrone blindspot as traced by empirical isochrones

Another interesting comparison that can be made using the empirical isochrones
of the archive concerns the isochrone blindspot between 100 and 500 Myr that
was initially introduced in Fig. 2. The isochrones corresponding to the lower
and upper age limits of said region are outlined in black in the panels of Fig
9 for the different _Gaia_ color combinations, again using the PARSEC code
with default parameters, no extinction and solar-like metallicity. Along with
the theoretical isochrones, the figure shows all empirical isochrones of the
archive clusters whose estimated ages coincide with the blindspot age range.

Figure 9: Visualization of empirical isochrones compared with theoretical ones
for the previously described isochrone blindspot ($\sim 100-500$ Myr). The
_black lines_ correspond to model isochrones delineating the edges of the
blindspot region, using PARSEC with default parameters. The _colored lines_
are all empirical isochrones belonging to clusters with estimated ages in this
age range.

It can be observed that the region where the theoretical isochrones overlap to
form the isochrone blindspot is not as narrow when empirical isochrones are
concerned. This broadening is most likely caused by the influence of the
extinction and metallicity values for each cluster, thereby again highlighting
the high dependence of the model-dependent isochrone fitting method on these
two input parameters. Another reason for the observed broadening could be the
uncertainty of the empirical isochrones around the main sequence turn-off
locations of the clusters at the upper end of the blindspot age range. In
cases where only a few select sources are observed in the red giant region,
the empirical isochrones may divert from the main sequence earlier than
theoretically predicted for the corresponding age. Finally, it is also
possible that some of the oldest clusters in the figure may have their ages
underestimated in literature, as their empirical main sequence turn-off points
seem lower than theory would suggest.

Nonetheless, we report a region between $\sim 5-9$ absolute G magnitudes,
where the empirical isochrones emulate their theoretical counterparts in
forming a very narrow section of overlapping curves – in a reduced form, the
blindspot remains visible, even using empirical isochrones.

Apart from the broadening of the upper main sequence, there is another
critical feature mapped by the empirical isochrones in the plot, which
concerns the lower main sequence (below the _dash-dotted line_): Both the
model isochrones and the empirical ones start to diverge from one another
again at this point. However, the empirical curves systematically assume a
different shape and curvature than the models, which is especially pronounced
for the _Gaia_ $\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ and
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}$ color indices. The color coding again
shows an evolutionary sequence for the empirical lines. The discrepancy
between models and observations is discussed further in Section 6.2.1.

###  5.3 Performance on different photometric systems

Figure 10: Empirical isochrones calculated for the clusters IC 4665 and the
Pleiades for six CMD variations using DANCe data.

To gauge the aspired flexibility and adaptive properties of the PCA-SVR
algorithm toward different CMD shapes and photometric systems other than
_Gaia_ , we applied our method to the DANCe dataset of the two clusters IC
4665 and the Pleiades (Olivares et al. 2018a; Miret-Roig et al. 2019). Figure
10 displays the results for six different passband combinations. The
isochrones calculated from these CMDs are also available online (see Appendix
B for details).

As can be seen, the algorithm demonstrated an impressive performance on all
tested CMDs, indicating great flexibility toward different kinds of input
data. Even for the mostly vertical CMD and characteristic knee-shaped kink of
the distribution of the $\mathrm{M}_{\mathrm{J}}$ vs. $\mathrm{J}-\mathrm{K}$
combination, reliable empirical isochrones that recover the distribution shape
in detail were found.

Compared to most of the results using _Gaia_ data, the edges of the
uncertainty bounds and the isochrones are considerably less well-behaved for
the new photometric systems, even describing abrupt turns in some cases.
However, they penetrate the very faint regions close to the observation limit,
which is why few sources are available for the fit. At the same time, the data
scatter is considerably more prominent than in most _Gaia_ CMDs. In light of
these complications, the performance of the algorithm is commendable. Even
more, an adaptation of the hyperparameters to the new photometric systems was
not necessary.

##  6 Discussion

In this section, we first address general quality aspects of using and
interpreting our empirical isochrones and possible limitations of the
algorithm before outlining concrete application scenarios of the archive as
age-scaling ladder or age brackets. As most of the new information gained from
empirical isochrones concerns the lower main sequence, we mainly focus on this
region. Furthermore, the performance and information gain of empirical
isochrones are compared with several popular models.

###  6.1 Empirical isochrone quality

The quality of the extracted empirical isochrones is evaluated in a two-
pronged manner: First, via comparison to the observed stellar distribution in
the CMD and considering the shape and width of the uncertainty bounds and
second, via the more general viewpoint regarding the sensitivity of the
extracted curves on variations in the key observational parameters from which
the CMD is created.

Regarding the former, it can be seen from the cluster CMDs in Figs. 15 and 16
that the isochrone extraction with PCA and SVR yields accurate representations
of the apparent shape of the cluster distributions. In general, there is
little fluctuation between the various resampled isochrones, and the
uncertainty regions indicated by _black lines_ are very narrow. The results
show a successful hyperparameter tuning and the general good applicability of
the machine learning procedure to the extraction task. The central parts of
the cluster main sequences are exceptionally well-defined for almost all
archive clusters. The only regions where the uncertainty bounds are larger are
the edge regions.

Regarding the latter, we identified five key parameters that may impact the
CMD of a stellar population itself, and by extension also the empirical
isochrone: the measurement uncertainties in photometry and parallaxes,
fraction of unresolved binaries, (differential) extinction, and field
contamination. We performed a detailed exploration of the sensitivity of our
proposed method to those parameters in Appendix E. To this end, we defined a
parameter range holding three values for each key parameter and constructed a
full factorial parameter grid consisting of 243 combinations. We then
evaluated each grid point on three representative archive clusters by creating
synthetic CMDs that incorporated the respective uncertainties in the key
parameters at the given grid point and extracting empirical isochrones from
them. We calculated an average deviation score between newly calculated and
original isochrones using a nearest-neighbor distance metric. By heuristically
determining thresholds for acceptable deviations between original and re-
computed isochrones, we determined parameter ranges for reliable results based
on our analysis (see Sect. E.2). We briefly summarize our findings in the
following paragraphs.

####  6.1.1 Edge regions

For a few clusters, the upper ends of the empirical isochrones start to
diverge and flatten toward the left of the populated area instead of matching
the brightest stars of the clusters. Prominent examples include the sigma Sco,
Antares, nu Cen, UPK 422, and NGC 2516 clusters in Figs. 15 and 16,
respectively. This usually happens for a specific combination of CMD
properties: Firstly, the upper main sequence only includes few and sparse
sources, creating an exceptionally high gradient in population density between
the upper and the lower mass end. Consequently, the bootstrapping process is
even less likely to draw sources representing the upper main sequence than the
general bias that the initial mass function imposes. Secondly, the brightest
sources are associated with larger errors, usually in their parallax
measurement, and therefore smaller weights in the regression. As a result, the
machine learning algorithm tends to avoid those data points, but as it is
designed to create a curve encompassing the entire dynamical range of the
cluster, it creates the previously described artifacts. Often, these are
already indicated by larger uncertainty bounds, but there are exceptions (see,
e.g., UPK 422). When working with such isochrones, one should consider
removing the empirical line’s final entry to have a more reliable
representation of the cluster.

####  6.1.2 Measurement uncertainties in photometry and parallax

Our sensitivity analysis showed, that neither photometric nor parallax
measurement uncertainties have significant influence on the extracted
isochrone. We purposefully used extreme assumptions for the errors in
attributing the chosen uncertainties to every star in a cluster instead of
drawing different values from a Gaussian distribution. Even so, we found that
viewed individually, their respective impact on the isochronal shape is
negligible. For the photometric uncertainty, we report an average deviation
between old and re-computed isochrone of less than 1 % for photometric
uncertainties $\leq 0.03$ mag for each star. The same holds true for a
fractional parallax error of 5 % for each star. For a fractional parallax
error of 10 %, the average deviation between the original and a re-computed
isochrone is around 2 % using a
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ passband combination in
the CMD, ca. 1 % using a $\mathrm{G}_{\mathrm{BP}}-\mathrm{G}$ CMD and below
one percent for the $\mathrm{G}-\mathrm{G}_{\mathrm{RP}}$ CMD (Tab. 13).
Nevertheless, we again note that in case of fractional parallax errors larger
than 10 %, our approach to distance estimation via parallax inversion would
not be applicable anymore.

When assuming the simultaneous presence of uncertainties in all parameters, we
could extract reliable isochrones for more than 50 % of the tested parameter
combinations for parallax uncertainties up to 5 % in the
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ color combination. For the
photometric uncertainties and other passband combinations, the fraction was
typically lower than 50 %.

####  6.1.3 Unresolved binary sequence and field star contamination

Concerning field star contaminants and the unresolved binary sequence, our
extraction method, in theory, should be robust against outliers as long as
their number does not exceed that of (resolved) cluster members. For the
archive clusters, this was only the case for the discarded RSG 7 cluster. We
achieve this robustness by regulating the $\varepsilon$-tube of the SVR
(Appendix D.2).

We verified this assumption in our sensitivity analysis in Appendix E by
including both the unresolved binary fraction and field star contamination as
parameters in our grid. When averaging over all three clusters and viewing the
parameter influence in isolation, meaning with all other parameter influences
set to zero (Tab. 13), we verified that 30 % binaries pose no problem to the
reliability of the empirical isochrones in any of the passband combinations.
On the other hand, the presence of 50 % unresolved binary contamination exceed
the reliability threshold, an expected outcome given the theoretical
limitations of our method. As for the field fraction, contamination up to 50 %
has no impact on the reliability of the extracted isochrones, when acting in
isolation.

We found that in the presence of other parameter uncertainties, disregarding
photometric uncertainties, empirical isochrones for the
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ passband combination could
be produced for more than 50 % of the grid points up to an unresolved binary
fraction of 35 to 40 % for the three test clusters. For the other passband
combinations the thresholds are not reached for all clusters, thus not
permitting a general statement. Interestingly, we find that artificially
adding field contamination to the data can statistically have beneficial
effects on intermediate age and old clusters, as they mostly add to the same
locus as the cluster (lower) main sequence stars. Due to its dependence on the
cluster age, the influence of field contamination on the isochrone reliability
in the presence of other uncertainty factors cannot be simply generalized.

####  6.1.4 Extinction

Observational data can be affected by interstellar extinction, which would
shift the position of an empirical isochrone toward lower magnitudes and
redder colors. Such mispositioning could potentially cause problems concerning
the application of empirical isochrones as age-scaling ladders. However,
extinction can be determined by using several sets of extra measurements, such
as stellar spectra, inference methods such as _Gaia_ Apsis (Fouesneau et al.
2023), or as a model parameter, for example via isochrone fitting. As a
result, similarly to cluster ages, one can find a vast array of different
extinction values for the archive clusters when consulting various works of
literature (Bossini et al. 2019; Cantat-Gaudin & Anders 2020; Dias et al.
2021). Adjusting the empirical isochrones for extinction would therefore
contradict the central focus of the isochrone archive of remaining model-
independent and free of any priors, as we would have to place preference on
one of the different approaches for extinction determination.

Still, the effects of the extinction parameter on the isochrone quality in our
specific use case need to be evaluated. As a first step, we calculated the
extinction vector using the coefficients determined by Wang & Chen (2019).
Since they only provide values for the second _Gaia_ data release, we use this
data as well, as it is enough for a qualitative picture of the situation. Fig.
11 shows a CMD of all sources included in our selections of Catalog I and
Catalog III, respectively.777As in Fig. 1, the vertical line of data points in
the CMD of Catalog I only appears in their DR2 selection and vanishes when
using DR3 data. The extinction vector in the _Gaia_ DR2 G band is displayed
alongside the data. It aligns well with the orientation of the bulk of the
main sequence locus of stars in both test cases, meaning that any extinction
present in the individual clusters would only result in a slight downward
shift along the main sequence. We further note that according to literature
(Bossini et al. 2019; Cantat-Gaudin & Anders 2020; Ratzenböck et al. 2023),
the highest estimated extinction of a cluster in the archive corresponds to
$A_{\mathrm{G}}=0.6312$ mag ($A_{V}=0.8$ mag) for UBC 17a.888Dias et al.
(2021) report extinction values $A_{V}>0.8$ mag for NGC 1662 ($A_{V}=1.214$
mag), Alessi 20 ($A_{V}=1.002$ mag) and Collinder 350 ($A_{V}=0.981$ mag).
However, these values are between a factor of 1.67 to 2.76 higher than those
reported by both Cantat-Gaudin & Anders (2020); Bossini et al. (2019). Any
shift would therefore only have a maximum length of around two-thirds of the
indicated extinction vector.

Figure 11: Scope and direction of the extinction vector, visualized on _Gaia_
DR2 data of the Catalog I and Catalog III selections using coefficients
published by Wang & Chen (2019). For both cases, the extinction vector is
almost parallel to the general orientation of the main sequence.

We validated our assumption that the extinction values for the clusters in the
archive are negligible in our sensitivity analysis in Appendix E. We created
synthetic CMDs for three representative archive clusters, assuming a flat
extinction level of either $A_{\mathrm{G}}=0.25$ mag or $A_{\mathrm{G}}=0.75$
mag, and calculated the deviation between the original and the newly computed
isochrone and compared the results to empirically derived threshold values for
each CMD passband combination and cluster, respectively. Viewed in isolation,
we also found that for all but the $\mathrm{G}-\mathrm{G}_{\mathrm{RP}}$ CMD,
an extinction value of $A_{\mathrm{G}}=0.25$ mag causes deviations well below
the threshold value; however, this does not hold true for the larger
extinction value of 0.75 mag (Tab. 13).

Given the presence of uncertainties in the other examined parameters, the
median isochronal deviation lies within the accepted range for
$A_{\mathrm{G}}\lesssim 0.3$ mag in the
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ CMD. For the other tow
combinations, only values well below 0.1 mag extinction should be considered
in the presence of further significant parametric uncertainties.

We also approximated the effects of differential extinction on the isochrone
extraction method (see Appendix E.3) by adding Gaussian noise to the flat
extinction values and fixing the other investigated parameters to plausible
values999Photometric uncertainty of 0.01 mag in all source, fractional
parallax uncertainty of 1% in all sources, 30 % unresolved binary fraction and
25 % field contamination fraction.. Given these assumptions, the reliability
threshold for a flat extinction assumption lies around 0.35 - 0.6 mag for the
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ CMD and around 0.2 - 0.3
mag for the $\mathrm{G}_{\mathrm{BP}}-\mathrm{G}$ CMD. We found that for
extinction levels $<0.5$ mag and standard deviations
$\sigma_{A_{\mathrm{G}}}\leq 1$ mag, the same thresholds apply as for the flat
extinction level.

We conclude, that for our archive clusters, we should be able to produce
reliable isochrones in the $\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$
passband combinations for almost all cases, even taking their estimated
extinction into account. The third quantile and maximum of our reference
extinction values are $A_{\mathrm{G}}=0.377$ mag and $A_{\mathrm{G}}=0.6312$
mag, respectively. The former lies within the reliability threshold range
defined in Appendix E.3, while the latter is slightly above the upper limit.
However, the actual uncertainties are likely to be below the assumptions of
the simulations, which are overestimating the errors especially for the
photometric and parallax uncertainties compared to a realistic distribution of
uncertainties. As the sensitivity to extinction seems to scale with cluster
age and/or distance, we searched the archive clusters for candidates with
similar properties as the test case NGC 752, meaning $d>400$ pc and an age of
$>1$ Gyr and checked whether their estimated extinction exceeds the
reliability threshold for NGC 752 of 0.3 mag. We then did the same search for
clusters similar to Blanco 1, that is $d>200pc$ and age $>100$ Myr, with an
extinction larger than the threshold value of 0.5 mag. None of the clusters in
the archive corresponding to said properties have reported extinction values
exceeding the respective threshold. The same applies to the
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}$ CMD thresholds. For the
$\mathrm{G}-\mathrm{G}_{\mathrm{RP}}$ passband combinations, we can only make
the definitive statement that we can produce reliable empirical isochrones for
flat extinction levels up to $A_{\mathrm{G}}=0.25$ mag, disregarding the
presence of uncertainties in other parameters.

It should be noted that users may subtract the extinction contribution from
the empirical isochrones themselves in the case of uniform extinction at any
time. In future works, we might also provide empirical isochrones considering
the extinction parameters inferred by _Gaia_ Apsis.

###  6.2 Lower main sequence analysis: Case studies

The lower main sequence of clusters is an important region for various
reasons. First and foremost, low-mass stars are known to be the most numerous
stars in the galaxy. Their properties have rendered them high-interest objects
in the search for exoplanets (see, e.g., Reiners et al. 2018), which in turn
makes estimations of their respective ages and masses paramount for furthering
our understanding, especially regarding planetary formation and habitability
research. Another vital aspect of why we should care for lower main sequence
stars is that it is often an essential discriminating feature of younger
populations ($\gtrsim 50-200$ Myr), as they often have no observable turn-off
yet, and fall into the blindspot region. As Fig. 9 showed, the lower main
sequence of such clusters can be empirically resolved and aid the correct age
ordering of such populations.

####  6.2.1 Theoretical vs. empirical isochrones down to the brown dwarf limit
for IC 4665

Modeling the dense and cool objects making up the faint end of (nearby)
stellar population CMDs with stellar evolution codes involves complex physics.
For instance, correct representations of nonideal effects within the stellar
equations of state, as well as for nuclear reaction processes, and accurate
stellar atmosphere models that account for the effects of molecular opacity
(Chabrier & Baraffe 1997; Baraffe et al. 1998) are needed. A particular
obstacle to the successful modeling of the evolution of low-mass objects is
posed by atmospheric convection or other stellar mixing processes such as
overshooting (Weiss & Heners 2013). Convection, at least in the mixing length
theory formalism, reportedly still has some limitations, for instance, in
mapping the thermal properties of convectively unstable atmospheres (Baraffe
et al. 2015). Despite the overall reported good agreement between models and
observations, especially in the optical bands, there is still often some
visible disagreement regarding the lower main sequence of stars that appears
over different models found in literature (e.g., Chen et al. 2014; Herczeg &
Hillenbrand 2015).

Figure 12: Comparison between theoretical and empirical isochrone performances
in the lower main sequence of the cluster IC 4665. _Left:_ Three theoretical
isochrones corresponding to an age of 30 Myr and using solar metallicity,
corrected regarding the estimated cluster extinction (Miret-Roig et al. 2019).
_Right:_ Empirical isochrones determined in this work (_red line_) and by
Miret-Roig et al. (2019) (_orange line_).

The disagreement between theory and observations for the faint end of the main
sequence has already been shown in Fig. 9 for the archive clusters. Now
instead, we want to demonstrate the sensitivity of stellar evolution models to
different input physics and how they compare to reality when applied
specifically to the very low-mass end of a population. To this end, we display
a $\mathrm{m}_{\mathrm{i}}$ vs. $\mathrm{i}-\mathrm{K}$ CMD of DANCe data of
the IC 4665 cluster in Fig. 12. We note that this plot displays apparent
magnitudes instead of the usual absolute magnitudes on its $y-$axis. The left
panel additionally displays model isochrones from BHAC15 (Baraffe et al.
2015), BT-Settl (Allard 2014), and PARSEC101010For this comparison, we chose
the OBC bolometric correction instead of the default choice of YBC+Vega, as
OBC matches the observations more closely. (Bressan et al. 2012),
corresponding to the age ($30$ Myr) and extinction ($A_{V}=0.72$ mag,
transformed into the used bands using tables from Wang & Chen 2019) values
used in Miret-Roig et al. (2019). It is immediately obvious that neither of
the three theoretical isochrones presents an accurate match for the data, but
even more astonishing is the fact that they are in such stark disagreement
with one another. The PARSEC isochrone does not penetrate the very faint
magnitudes and has the least resemblance to the observations. The BHAC15
isochrone does not match the population very well, either. The best fit is
achieved with the BT-Settl isochrone, but it still misses the sources at the
faint end. Our findings support the age variations between different models
already reported in the introduction in Sect. 2.1. But even more, we emphasize
that these isochrones are not only supposed to represent the same age but also
identical metallicity and extinction values.

In the right panel of the figure, the empirical isochrone derived by Miret-
Roig et al. (2019), and the one calculated in this work are displayed. Both
are much better fits than their theoretical counterparts for IC 4665. The
differences between our and the DANCe isochrones are a consequence of the
different extraction methods and the fact that the DANCe isochrone traces the
lower left edge of the distribution, whereas we decided to derive the curve
from the densest part, but the differences are less pronounced than for the
model isochrones. In the areas with less scatter in the data, both empirical
curves are in excellent agreement and accurately trace the sequence down to
the faintest magnitudes.

The case study on IC 4665 shows that our observational data are good enough to
map the observed distribution in detail, meaning that empirical isochrones
could be important for constraining stellar evolution models. We anticipate
that our results will provide reliable calibration data, specifically in the
lower mass regime, where the physics gets increasingly complex, albeit they
could also serve for calibration purposes along the rest of the main sequence.

####  6.2.2 Age-bracketing of the Meingast 1 stream

Figure 13: Age bracketing of the Meingast 1 stream with empirical isochrones,
whose corresponding cluster CMDs are shown in the right column of the figure.
The gray points indicate the combined data of all three clusters, and the
population connected to the plotted isochrone is overlaid in blue. By
comparing the empirical isochrone of the stellar stream (_yellow line_) with
similar archive isochrones, we bracketed the stream and, in extension, its
relative age, by the Pleiades (_magenta line_) and by the Blanco 1 cluster
(_green line_).

Figure 14: Comparison of the difference between the empirical isochrones in
the lower main sequence region with theoretical isochrones across a realistic
extinction and metallicity range._Left:_ PARSEC model isochrones with
different possible extinction and metallicity values found for the Pleiades
and the stream in literature. _Right:_ Zoom-in of the lower main sequence
region of the age bracket of the Meingast 1 stream shown in Fig. 13.

We illustrate the idea behind the relative age determination using empirical
isochrones on a type of stellar conglomerate discovered fairly recently – the
stellar stream Meingast 1 (Meingast et al. 2019). As a relatively unfamiliar
object with a well-defined CMD member selection and an estimated age right
inside the blindspot region, it presents a particularly interesting case study
for assessing the usefulness of the empirical isochrone archive.

The first discovered example of its class, the stream is a prime example of
how easily isochronal ages can be distorted by nonmember stars at crucial
points such as the turn-off region: Upon its discovery, the stream was
estimated to be around 1 Gyr (Meingast et al. 2019), due to the presence of a
single evolved star right at the top of the member selection in the CMD.
Subsequent works and age determinations using gyrochronology (Curtis et al.
2019), a lithium-rotation connection in low-mass members (Arancibia-Silva et
al. 2020), isochrone fitting with an updated member selection (Ratzenböck et
al. 2020; Röser & Schilbach 2020), lithium abundance measurements (Hawkins et
al. 2020), and empirical lithium equivalent widths (Jeffries et al. 2023) have
revised the original estimate and settled on an age around $110-120$ Myr.

However, apart from the discovery paper, all literature works have compared
the stream to the Pleiades cluster, and almost all estimates even rest almost
solely on this fact. For instance, its gyrochronological age was determined
based on its similarity to one of only three benchmark clusters referenced as
120 Myr, 670 Myr, and 1 Gyr, respectively. Similarly, the age estimate from
lithium abundance measurements for both the Pleiades and the stream has only
very broad constraints of $\sim 120-200$ Myr and stems from comparing the
stream data with measurements from the Pleiades, the Hyades, and the galactic
disk. While the results indeed show that an apparent similarity exists between
the two populations that can be traced by different indicators, a comparison
between so few objects, though necessitated by lack of further observations,
does not permit the determination of precise (relative) age limits.

Once again, one has to fall back onto model-dependent approaches to get more
precise age estimates. But for the Meingast 1 stream, even this strategy is
hampered as it, along with the Pleiades, lies squarely inside the isochrone
blindspot (Sect. 2.2). This is also reflected in literature, where it is
reported that PARSEC isochrones for the stream and the Pleiades show no
significant differences between the 80 and 130 Myr (Curtis et al. 2019), or
100 and 150 Myr (Röser & Schilbach 2020).

We, on the other hand, now have ample empirical information in the blindspot
region, specifically in the lower main sequence, thanks to the empirical
isochrone archive. Hence, we could select the ones best resembling the stream
and determine relative “age brackets” for the stream, which are shown in Fig.
13: The upper bracket is provided by the Pleiades, whereas the Blanco 1
cluster forms the lower bracket. For our subsequent analysis, we focus
exclusively on the regions where the uncertainty regions (_dashed lines_) and
the isochrones themselves are indistinguishable to minimize uncertainties
brought on by the quality limitations of the empirical isochrones (Sect. 6.1).

Similarly to the other literature works, we find a close resemblance between
the Pleiades and the stream over almost the entire blindspot region. But we
also see a clear shift of the empirical lower main sequences of the two
populations starting around $\mathrm{M}_{\mathrm{G}}\approx 9$ mag and
continuing until the broadening of the uncertainty bounds. Judging from only
the lower main sequence, the Pleiades would therefore be slightly younger than
Meingast 1. However, this is partly contested by literature (Table 2) and also
not supported by the upper main sequence of the empirical isochrones.

However, we find that along the lower main sequence the Meingast 1 stream
isochrone matches the one of the Blanco 1 cluster almost exactly. The stream
and Blanco 1 are both located below the galactic plane in the direction of the
galactic south pole, indicating that their extinction should be negligible. In
contrast, the Pleiades are often attributed a low but not zero extinction,
although its value, when derived via isochrone fitting, is affected by the
same uncertainties as the age estimates due to the strong parameter
degeneracy.

Another factor to consider in this difference between the isochrones is
metallicity: Investigations of Ratzenböck et al. (2020) indicate a slight
difference between the stream’s and the Pleiades’ metallicity values of
$\Delta Z\simeq 0.002$. To determine whether the shift between the empirical
isochrones of the two populations could be a result of either of the two
parameters, we calculate PARSEC isochrones using the highest estimated
extinction for the Pleiades found in our literature sources (Dias et al. 2021,
$A_{V}=0.18$ mag), and stellar metallicity values ranging from solar to
$Z=0.018$. The results of this qualitative analysis are displayed in Fig. 14.
As can be seen from the comparison between the left and the right panel, the
theoretical isochrones do not disperse enough to account for the apparent
differences in the isochrones. In other words, the shift in the lower main
sequence between the Pleiades and the stream cannot be explained solely by the
combined effects of extinction and metallicity.

Judging from the results of the empirical lower main sequence analysis, we
propose that the age of the Meingast 1 stream is between that of Blanco 1 and
the Pleiades, with slightly more resemblance to the former population, than to
the latter. Alas, comparing literature ages of those two populations yields a
disagreement on which one is actually older and which is younger (see Table
2). Again, this highlights how flimsy model age determination can be.

Table 2: Different age estimates for the clusters bracketing the Meingast 1 stream. Since there is a general disagreement on which of the clusters is the older one, providing a relative age estimate for the stream is the best possibility. Reference | Blanco 1 | Pleiades
---|---|---
_Isochronal log(age)_
Bossini et al. (2019) | 7.975 | 7.973
Cantat-Gaudin et al. (2020) | 8.02 | 7.89
Dias et al. (2021) | 8.012 | 8.116
_Empirical lithium EW log(age)_
Jeffries et al. (2023) | 7.86 | 7.97

Despite the disagreements in literature on which of the bracketing clusters
actually corresponds to the lower and upper age limit, respectively, we can
still provide a quite narrow relative age range of the stream. Given its
position relative to the Pleiades in the lower main sequence, it seems as if
it were a few Myr younger, even after factoring in possible extinction. From
this perspective, it also seems as if the Blanco 1 cluster were indeed older
than the Pleiades. This trend does seem reversed in the upper main sequence,
but we again highlight that the regions with diverging uncertainty bounds
should be interpreted with much more care than those with narrow uncertainty
regions.

##  7 Summary and conclusions

The new era of high-precision astrometric and photometric data provided in
large quantities by all-sky surveys such as _Gaia_ allows us to find new
synergies between astronomy and data science. In doing so, we were able to add
new, purely observation-based information to the complex topic of stellar age
determination.

We developed a fast, flexible algorithm based on machine learning and
statistical analysis tools that produces empirical isochrones for a multitude
of CMD shapes and evolutionary stages of clusters. Its only inputs are the
photometric and astrometric data needed to create any arbitrary color-
magnitude diagram of a cluster or stellar population. By applying our
algorithm to different recently published cluster member catalogs, we created
an archive of empirical isochrones that includes 83 clusters within 500 pc
representing an age range between 7 Myr and 3 Gyr. We draw the following
conclusions from our results and subsequent analyses:

  1. 1.

When observing the general, larger trends, isochronal age determination is in
good agreement with the empirical isochrone positions (Fig. 7).

  2. 2.

Validation against the Sco-Cen data of Catalog III (Fig. 8) shows that it is
possible to create a fine grid for relative age determinations and that for
homogeneously determined groups and ages there exists a one-to-one correlation
between the estimated absolute ages and the distribution of the empirical
isochrones.

  3. 3.

The algorithm works for different photometric systems and CMD configurations
(Fig. 10), even though it was developed using only _Gaia_ data.

  4. 4.

The empirical isochrones are overall detailed and accurate CMD
representations. High data scatter and sudden source drops are the limiting
factors in the isochrone quality (Sect. 6.1).

  5. 5.

A detailed sensitivity analysis concluded that, viewed in isolation, the
empirical isochrones are reliable against at least a photometric measurement
uncertainties of 0.03 mag, a fractional parallax uncertainties of 10 %, a
binary contamination of 30 %, an extinction level of 0.25 mag
$A_{\mathrm{G}}$111111Except for the $\mathrm{G}-\mathrm{G}_{\mathrm{RP}}$
passband combination., and a field contamination fraction of 50 %. A more
detailed analysis of the effects of extinction concluded, that under the
assumption of 0.01 mag photometric uncertainty, 1 % parallax uncertainty, 30 %
unresolved binary fraction and 25 % field contamination, empirical isochrones
remain reliable given a mean extinction of 0.3 - 0.6 mag, depending on the
cluster age, in the $\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ CMD,
depending on the inspected cluster. For the
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}$ CMD, this value range shifts to 0.2 -
0.3 mag. When considering combined uncertainty effects, the reliability may
vary depending on the specific uncertainty values (Sect. E).

  6. 6.

Empirical isochrones add considerable information toward the hard-to-model
lower main sequence. This is especially important in the isochrone blindspot
region.

  7. 7.

There is still ample room for improvement when comparing model isochrones with
observations, specifically at low stellar masses (Sect. 6.2). We are aware
that modeling the physics of low-mass stars is particularly challenging and
anticipate that our empirical isochrones could serve for calibration purposes
in this regard (Fig. 12).

  8. 8.

The archive of isochrones can be used as age-scaling ladder for determining
relative ages via empirical isochrone bracketing. We showed this in a case
study of the newly found Meingast 1 stream, which is closely bracketed by not
only the Pleiades cluster as established in the literature but also shares a
lot of resemblance with the Blanco 1 cluster (Fig. 13).

In the future, we aim to expand the empirical isochrone archive toward an even
larger sample of nearby open clusters and look for loci of empirical lines
that can be connected to similar global attributes of open clusters.

###### Acknowledgements.

The authors thank the referee Javier Olivares for providing very useful
comments that helped to improve the quality of this work.The authors want to
thank the research group in Vienna and especially João Alves, Josefa
Großschedl, Núria Miret-Roig and Sebastian Ratzenböck, for their valuable
insights, input and discussion to this work. Funded by the European Union
(ERC, ISM-FLOW, 101055318). Views and opinions expressed are however those of
the author(s) only and do not necessarily reflect those of the European Union
or the European Research Council Executive Agency. Neither the European Union
nor the granting authority can be held responsible for them. 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. This research made use of
Astropyhttp://www.astropy.org, a community-developed core Python package for
Astronomy (Astropy Collaboration et al. 2013, 2018). We also acknowledge the
use of Python https://www.python.org, along with the packages that were used
in the data analysis of this work, including NumPy (van der Walt et al. 2011),
scikit-learn (Pedregosa et al. 2011), Matplotlib (Hunter 2007), and Plotly
(Plotly Technologies Inc. 2015). This research has made use of the VizieR
catalog access tool, CDS, Strasbourg, France (Ochsenbein et al. 2000; DOI:
10.26093/cds/vizier). This research has made use of TOPCAT, an interactive
graphical viewer and editor for tabular data (Taylor 2005).

## References

  * Adamo et al. (2013) Adamo, A., Östlin, G., Bastian, N., et al. 2013, ApJ, 766, 105
  * Allard (2014) Allard, F. 2014, in Exploring the Formation and Evolution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham, Vol. 299, 271–272
  * Arancibia-Silva et al. (2020) Arancibia-Silva, J., Bouvier, J., Bayo, A., et al. 2020, A&A, 635, L13
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, Astronomical Journal, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, Astronomy and Astrophysics, 558, A33
  * Bailer-Jones (2015) Bailer-Jones, C. A. L. 2015, PASP, 127, 994
  * Baraffe et al. (1998) Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42
  * Barnes (2007) Barnes, S. A. 2007, ApJ, 669, 1167
  * Basri et al. (1996) Basri, G., Marcy, G. W., & Graham, J. R. 1996, ApJ, 458, 600
  * Bell et al. (2012) Bell, C. P. M., Naylor, T., Mayne, N. J., Jeffries, R. D., & Littlefair, S. P. 2012, MNRAS, 424, 3178
  * Binks & Jeffries (2014) Binks, A. S. & Jeffries, R. D. 2014, MNRAS, 438, L11
  * Bishop (2006) Bishop, C. M. 2006, Pattern Recognition and Machine Learning (Information Science and Statistics) (Berlin, Heidelberg: Springer-Verlag)
  * Bossini et al. (2019) Bossini, D., Vallenari, A., Bragaglia, A., et al. 2019, Astronomy and Astrophysics, 623, A108
  * Bouy et al. (2013) Bouy, H., Bertin, E., Moraux, E., et al. 2013, A&A, 554, A101
  * Bouy et al. (2015) Bouy, H., Bertin, E., Sarro, L. M., et al. 2015, Astronomy and Astrophysics, 577, A148
  * Brandner et al. (2023) Brandner, W., Calissendorff, P., & Kopytova, T. 2023, arXiv e-prints, arXiv:2306.03132
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, Monthly Notices of the Royal Astronomical Society, 427, 127
  * Cantat-Gaudin & Anders (2020) Cantat-Gaudin, T. & Anders, F. 2020, Astronomy and Astrophysics, 633, A99
  * Cantat-Gaudin et al. (2020) Cantat-Gaudin, T., Anders, F., Castro-Ginard, A., et al. 2020, A&A, 640, A1
  * Cantat-Gaudin et al. (2018) Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018, Astronomy and Astrophysics, 618, A93
  * Castro-Ginard et al. (2022) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2022, A&A, 661, A118
  * Castro-Ginard et al. (2018) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2018, A&A, 618, A59
  * Castro-Ginard et al. (2021) Castro-Ginard, A., McMillan, P. J., Luri, X., et al. 2021, A&A, 652, A162
  * Chabrier & Baraffe (1997) Chabrier, G. & Baraffe, I. 1997, A&A, 327, 1039
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560
  * Chen et al. (2015) Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068
  * Chen et al. (2014) Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525
  * Chevance et al. (2020) Chevance, M., Kruijssen, J. M. D., Vazquez-Semadeni, E., et al. 2020, Space Sci. Rev., 216, 50
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Cortes & Vapnik (1995) Cortes, C. & Vapnik, V. 1995, Machine Learning, 20, 273
  * Curtis et al. (2019) Curtis, J. L., Agüeros, M. A., Mamajek, E. E., Wright, J. T., & Cummings, J. D. 2019, AJ, 158, 77
  * Da Rio et al. (2017) Da Rio, N., Tan, J. C., Covey, K. R., et al. 2017, Astrophysical Journal, 845, 105
  * David & Hillenbrand (2015) David, T. J. & Hillenbrand, L. A. 2015, Astrophysical Journal, 804, 146
  * de Boer & Seggewiss (2008) de Boer, K. & Seggewiss, W. 2008, Stars and Stellar Evolution (Les Ulis: EDP Sciences)
  * De Rosa et al. (2014) De Rosa, R. J., Patience, J., Wilson, P. A., et al. 2014, MNRAS, 437, 1216
  * Dias et al. (2021) Dias, W. S., Monteiro, H., Moitinho, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 504, 356
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8
  * Draine (2003) Draine, B. T. 2003, ApJ, 598, 1017
  * Ester et al. (1996) Ester, M., Kriegel, H.-P., Sander, J., & Xu, X. 1996, in Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, KDD’96 (AAAI Press), 226–231
  * Feiden (2016) Feiden, G. A. 2016, A&A, 593, A99
  * Fouesneau et al. (2023) Fouesneau, M., Frémat, Y., Andrae, R., et al. 2023, A&A, 674, A28
  * Fürnkranz et al. (2019) Fürnkranz, V., Meingast, S., & Alves, J. 2019, A&A, 624, L11
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018a, Astronomy and Astrophysics, 616, A10
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018b, Astronomy and Astrophysics, 616, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, Astronomy and Astrophysics, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, Astronomy and Astrophysics, 595, A1
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211
  * Gaia Collaboration et al. (2017) Gaia Collaboration, van Leeuwen, F., Vallenari, A., et al. 2017, A&A, 601, A19
  * Hastie & Stuetzle (1989) Hastie, T. & Stuetzle, W. 1989, Journal of the American Statistical Association, 84, 502
  * Hawkins et al. (2020) Hawkins, K., Lucey, M., & Curtis, J. 2020, MNRAS, 496, 2422
  * Herczeg & Hillenbrand (2015) Herczeg, G. J. & Hillenbrand, L. A. 2015, ApJ, 808, 23
  * Hidalgo et al. (2018) Hidalgo, S. L., Pietrinferni, A., Cassisi, S., et al. 2018, ApJ, 856, 125
  * Hosek et al. (2020) Hosek, Matthew W., J., Lu, J. R., Lam, C. Y., et al. 2020, AJ, 160, 143
  * Hotelling (1933) Hotelling, H. 1933, Journal of educational psychology, 24, 417
  * Hunt & Reffert (2023) Hunt, E. L. & Reffert, S. 2023, arXiv e-prints, arXiv:2303.13424
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Hurley & Tout (1998) Hurley, J. & Tout, C. A. 1998, Monthly Notices of the Royal Astronomical Society, 300, 977
  * Jeffries et al. (2023) Jeffries, R. D., Jackson, R. J., Wright, N. J., et al. 2023, MNRAS, 523, 802
  * Jolliffe (2002) Jolliffe, I. T. 2002, Principal component analysis for special types of data (Springer)
  * Jolliffe & Cadima (2016) Jolliffe, I. T. & Cadima, J. 2016, Philosophical Transactions of the Royal Society of London Series A, 374, 20150202
  * Kerr et al. (2022a) Kerr, R., Kraus, A. L., Murphy, S. J., et al. 2022a, ApJ, 941, 143
  * Kerr et al. (2022b) Kerr, R., Kraus, A. L., Murphy, S. J., et al. 2022b, ApJ, 941, 49
  * Kerr et al. (2021) Kerr, R. M. P., Rizzuto, A. C., Kraus, A. L., & Offner, S. S. R. 2021, ApJ, 917, 23
  * Krone-Martins & Moitinho (2014) Krone-Martins, A. & Moitinho, A. 2014, Astronomy and Astrophysics, 561, A57
  * Lada & Lada (2003) Lada, C. J. & Lada, E. A. 2003, Annual Reviews of Astronomy and Astrophysics, 41, 57
  * Li & Shao (2022) Li, L. & Shao, Z. 2022, ApJ, 930, 44
  * Li et al. (2020) Li, L., Shao, Z., Li, Z.-Z., et al. 2020, ApJ, 901, 49
  * Li et al. (2021) Li, Z.-Z., Li, L., & Shao, Z. 2021, Astronomy and Computing, 36, 100483
  * Lindegren et al. (2018) Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2
  * Mamajek & Bell (2014) Mamajek, E. E. & Bell, C. P. M. 2014, MNRAS, 445, 2169
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77
  * Martín et al. (2018) Martín, E. L., Lodieu, N., Pavlenko, Y., & Béjar, V. J. S. 2018, ApJ, 856, 40
  * Meingast et al. (2019) Meingast, S., Alves, J., & Fürnkranz, V. 2019, A&A, 622, L13
  * Meingast et al. (2021) Meingast, S., Alves, J., & Rottensteiner, A. 2021, Astronomy and Astrophysics, 645, A84
  * Messina et al. (2022) Messina, S., Nardiello, D., Desidera, S., et al. 2022, A&A, 657, L3
  * Miret-Roig et al. (2024) Miret-Roig, N., Alves, J., Barrado, D., et al. 2024, Nature Astronomy, 8, 216
  * Miret-Roig et al. (2019) Miret-Roig, N., Bouy, H., Olivares, J., et al. 2019, Astronomy and Astrophysics, 631, A57
  * Miret-Roig et al. (2020) Miret-Roig, N., Galli, P. A. B., Brandner, W., et al. 2020, A&A, 642, A179
  * Murray (2011) Murray, N. 2011, ApJ, 729, 133
  * Netopil et al. (2015) Netopil, M., Paunzen, E., & Carraro, G. 2015, A&A, 582, A19
  * Nguyen et al. (2022) Nguyen, C. T., Costa, G., Girardi, L., et al. 2022, A&A, 665, A126
  * Ochsenbein et al. (2000) Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23
  * Olivares et al. (2023) Olivares, J., Bouy, H., Miret-Roig, N., et al. 2023, A&A, 671, A1
  * Olivares et al. (2019) Olivares, J., Bouy, H., Sarro, L. M., et al. 2019, Astronomy and Astrophysics, 625, A115
  * Olivares et al. (2021) Olivares, J., Bouy, H., Sarro, L. M., et al. 2021, A&A, 649, A159
  * Olivares et al. (2018a) Olivares, J., Sarro, L. M., Moraux, E., et al. 2018a, A&A, 617, A15
  * Olivares et al. (2018b) Olivares, J., Sarro, L. M., Moraux, E., et al. 2018b, A&A, 617, A15
  * Pamos Ortega et al. (2022) Pamos Ortega, D., García Hernández, A., Suárez, J. C., et al. 2022, MNRAS, 513, 374
  * Pearson (1901) Pearson, K. F. 1901, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2, 559
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Pelkonen et al. (2024) Pelkonen, V. M., Miret-Roig, N., & Padoan, P. 2024, A&A, 683, A165
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49
  * Pietrinferni et al. (2021) Pietrinferni, A., Hidalgo, S., Cassisi, S., et al. 2021, ApJ, 908, 102
  * Plotly Technologies Inc. (2015) Plotly Technologies Inc. 2015, Collaborative data science
  * Ratzenböck et al. (2023) Ratzenböck, S., Alves, J., Großschedl, J. E., et al. 2023, arXiv e-prints, arXiv:2302.07853
  * Ratzenböck et al. (2022) Ratzenböck, S., Großschedl, J. E., Möller, T., et al. 2022, arXiv e-prints, arXiv:2211.14225
  * Ratzenböck et al. (2020) Ratzenböck, S., Meingast, S., Alves, J., Möller, T., & Bomze, I. 2020, A&A, 639, A64
  * Reiners et al. (2018) Reiners, A., Zechmeister, M., Caballero, J. A., et al. 2018, A&A, 612, A49
  * Röser & Schilbach (2020) Röser, S. & Schilbach, E. 2020, A&A, 638, A9
  * Salaris et al. (2022) Salaris, M., Cassisi, S., Pietrinferni, A., & Hidalgo, S. 2022, MNRAS, 509, 5197
  * Sarro et al. (2014) Sarro, L. M., Bouy, H., Berihuete, A., et al. 2014, A&A, 563, A45
  * Sim et al. (2019) Sim, G., Lee, S. H., Ann, H. B., & Kim, S. 2019, Journal of Korean Astronomical Society, 52, 145
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Smola & Schölkopf (2004) Smola, A. J. & Schölkopf, B. 2004, Statistics and Computing, 14, 199–222
  * Soderblom (2010) Soderblom, D. R. 2010, Annual Reviews of Astronomy and Astrophysics, 48, 581
  * Soderblom et al. (2014) Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E., & Naylor, T. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 219
  * Squicciarini & Bonavita (2022) Squicciarini, V. & Bonavita, M. 2022, A&A, 666, A15
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Tipping & Bishop (1999) Tipping, M. E. & Bishop, C. M. 1999, Neural Computation, 11, 443
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engg., 13, 22
  * Wang & Chen (2019) Wang, S. & Chen, X. 2019, ApJ, 877, 116
  * Weiss & Heners (2013) Weiss, A. & Heners, N. 2013, in European Physical Journal Web of Conferences, Vol. 43, European Physical Journal Web of Conferences, 01002

##  Appendix A Cluster CMDs

Figure 15: CMDs of the individual archive clusters (1-40), along with their
calculated empirical isochrones and uncertainty bounds. The clusters are
ordered according to their literature ages in ascending fashion.

Figure 16: Continuation of Fig. 15 for the remaining archive clusters (41-83).

To demonstrate the general quality of the empirical isochrones derived in this
work, we showcase the results for all 83 clusters in our selection in the
_Gaia_ $\mathrm{G}$ vs. $\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$
CMDs in Figs. 15 and 16. The clusters are arranged in ascending order based on
their estimated ages (see Table 1). The empirical isochrone is shown in red,
while the 5th and 95th percentile boundaries, calculated from 1000
bootstrapped curves, are represented by _black lines_. It can be seen that
while there is generally a good agreement between the literature ages and the
evolution of the empirical isochrone shapes, there are also some obvious
differences in the ordering, especially when inspecting rather evolved
clusters with stars populating the giant region.

##  Appendix B Empirical isochrone tables

In the following, we show an excerpt of one of the result tables for the
empirical isochrones, as it may be found online. Table 12 displays the first
ten entries of the first archive cluster for the _Gaia_ DR3
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ passbands (the full table
is available online). The same information is available for the
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ and
$\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$ passband variations for
all archive clusters. There are two distinct isochrones per passband
combination available for the ten clusters appearing in both Catalog I and
Catalog II, between we distinguish by assigning the empirical isochrones
corresponding to the second catalog the suffix “C2.” Our adopted reference age
estimates in this manuscript can be found in the last column of Table 12.
Similar tables are available for the same passband combinations using _Gaia_
DR2 and for the DANCe CMDs shown in Fig. 10, respectively.

Table 3: First ten table entries sorted by cluster name, showing only the
columns concerning the $\mathrm{G}_{\mathrm{BP}}-\mathrm{G}_{\mathrm{RP}}$
combination of _Gaia_ DR3, as well as the reference age used in this work. The
columns correspond to the lower (lb) and upper (ub) uncertainty bounds, as
well as the empirical isochrone (iso) calculated for each cluster with
$n_{\mathrm{boot}}=1000$.

Cluster | BPRP_lb_x | BPRP_lb_y | BPRP_iso_x | BPRP_iso_y | BPRP_ub_x | BPRP_ub_y | … | ref_age
---|---|---|---|---|---|---|---|---
ASCC_127 | 0.0311 | 0.5926 | -0.0110 | 0.6027 | -0.3193 | 0.6772 | … | 7.26
ASCC_127 | 0.0876 | 1.2667 | 0.0397 | 1.2782 | -0.0328 | 1.2957 | … | 7.26
ASCC_127 | 0.1387 | 1.5687 | 0.0660 | 1.5863 | 0.0296 | 1.5951 | … | 7.26
ASCC_127 | 0.2099 | 2.1380 | 0.1213 | 2.1594 | 0.0958 | 2.1655 | … | 7.26
ASCC_127 | 0.2123 | 2.1846 | 0.1260 | 2.2055 | 0.1004 | 2.2116 | … | 7.26
ASCC_127 | 0.2165 | 2.2862 | 0.1362 | 2.3056 | 0.1113 | 2.3116 | … | 7.26
ASCC_127 | 0.2472 | 2.8702 | 0.2128 | 2.8785 | 0.1932 | 2.8832 | … | 7.26
ASCC_127 | 0.2517 | 2.9404 | 0.2226 | 2.9475 | 0.2038 | 2.9520 | … | 7.26
ASCC_127 | 0.2610 | 3.0639 | 0.2366 | 3.0698 | 0.2196 | 3.0739 | … | 7.26
ASCC_127 | 0.2767 | 3.1764 | 0.2504 | 3.1828 | 0.2318 | 3.1873 | … | 7.26
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | … | ⋮

121212The full table will be made available at the CDS.

##  Appendix C Preprocessing the data with PCA

Although the concept of PCA was already developed more than a century ago,
(Pearson 1901; Hotelling 1933), its practical feasibility only started
thriving with the emergence of computers. Since then, PCA has received
recognition primarily as a dimensionality reduction method that unifies the
ideas of increasing data interpretability and minimal information loss. Its