# Exploring the Chemistry and Mass Function of the Globular Cluster 47 Tucanae
with New Theoretical Color-Magnitude Diagrams

Roman Gerasimov Center for Astrophysics and Space Sciences, University of
California, San Diego, La Jolla, California 92093, USA  Department of Physics
and Astronomy, University of Notre Dame, Nieuwland Science Hall, Notre Dame,
IN 46556, USA  Adam J. Burgasser Center for Astrophysics and Space Sciences,
University of California, San Diego, La Jolla, California 92093, USA  Ilaria
Caiazzo Division of Physics, Mathematics and Astronomy, California Institute
of Technology, Pasadena, CA 91125, USA  Derek Homeier Aperio Software Ltd.,
Insight House, Riverside Business Park, Stoney Common Road, Stansted, Essex,
CM24 8PL, United Kingdom  Harvey B. Richer Department of Physics and
Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver,
BC V6T 1Z1, Canada  Matteo Correnti INAF Osservatorio Astronomico di Roma, Via
Frascati 33, 00078, Monteporzio Catone, Rome, Italy  ASI-Space Science Data
Center, Via del Politecnico, I-00133, Rome, Italy  Jeremy Heyl Department of
Physics and Astronomy, University of British Columbia, 6224 Agricultural Road,
Vancouver, BC V6T 1Z1, Canada

###### Abstract

Despite their shared origin, members of globular clusters display star-to-star
variations in composition. The observed pattern of element abundances is
unique to these stellar environments, and cannot be fully explained by any
proposed mechanism. It remains unclear whether stars form with chemical
heterogeneity, or inherit it from interactions with other members. These
scenarios may be differentiated by the dependence of chemical spread on
stellar mass; however, obtaining a sufficiently large mass baseline requires
abundance measurements on the lower main sequence that is too faint for
spectroscopy even in the nearest globular clusters. We developed a stellar
modelling method to obtain precise chemical abundances for stars near the end
of the main sequence from multiband photometry, and applied it to the globular
cluster 47 Tucanae. The computational efficiency is attained by matching
chemical elements to the model components that are most sensitive to their
abundance. We determined [O/Fe] for $\sim 5000$ members below the main
sequence knee at the level of accuracy, comparable to the spectroscopic
measurements of evolved members in literature. The inferred distribution
disfavors stellar interactions as the origin of chemical spread; however, an
accurate theory of accretion is required to draw a more definitive conclusion.
We anticipate that future observations of 47 Tucanae with JWST will extend the
mass baseline of our analysis into the substellar regime. Therefore, we
present predicted color-magnitude diagrams and mass-magnitude relations for
the brown dwarf members of 47 Tucanae.

Globular clusters (656) — Stellar populations (1622) — Chemical abundances
(224) — Brown dwarfs (185) — HST photometry (756)

tablenum SIXtablenum

##  1 Introduction

Globular clusters are large gravitationally bound conglomerations of stars,
commonly observed in the stellar halos of galaxies including our own (Harris,
1991; Harris & Racine, 1979; Hubble, 1932). Since the brightest examples (e.g.
$\omega$ Centauri, 47 Tucanae) are prominent naked-eye targets, the study of
globular clusters dates back multiple centuries (de Lacaille, 1755; Messier,
1781). A typical globular cluster hosts over ${10}^{5}$ members (Hilker et
al., 2020), far outnumbering the populations of most open clusters ($\sim 20$
– ${10}^{4}$ members, Qin et al. 2023; Sampedro et al. 2017) that form in the
present-day Milky Way. This points to the origin of globular clusters in the
early phases of the hierarchical assembly of galaxies (Freeman & Bland-
Hawthorn, 2002), making them some of the oldest objects in the universe, often
exceeding $\qty{12}{Gyr}$ in age (Carretta et al., 2000; VandenBerg et al.,
2013). By consequence, the members of most globular clusters exhibit sub-solar
metallicities ($\ce{[Fe/H]}\lesssim\qty{-1}{dex}$, Bailin 2019; Harris 1996;
Baum 1952) and have a distinct color distribution (ten Bruggencate, 1927;
Hachenberg, 1939; Greenstein, 1939) that motivated the original definition of
stellar populations (Baade, 1944; Oort, 1926).

Globular clusters are prime stellar astrophysics laboratories (Moehler, 2001)
and probes of galactic formation and evolution (Lee, 2012; Recio-Blanco, 2018;
Arkelyan & Pilipenko, 2022; Strader et al., 2004). The old ages of globular
clusters provide a direct constraint on the age of the universe (Krauss &
Chaboyer, 2003; Valcin et al., 2020, 2021), while their kinematic properties
trace the distribution of dark matter (Doppel et al., 2021; Alabi et al.,
2016, 2017). Furthermore, tidally disrupted globular clusters are responsible
for the production of stellar streams, alongside dwarf galaxies (Newberg,
2016). Understanding the properties and evolution of these unique objects is
therefore of uttermost importance to stellar astrophysics, galactic science
and cosmology.

Naively, the coeval nature of star clusters leads to the expectation of
uniformity in element abundances among their members. Nonetheless, the
presence of star-to-star chemical variations within globular clusters was
evident since the early spectroscopic observations that found differing
cyanogen (CN) absorption strength (Popper, 1947) in stars of similar spectral
types and luminosity classes (the “cyanogen discrepancy”, Keenan & Keller
1953), followed by the first discoveries of carbon star members (Harding,
1962; Wing & Stock, 1973). At first, these anomalies were largely attributed
to evolutionary effects rather than a genuine primordial inhomogeneity (e.g.,
Zinn 1973; Norris & Bessell 1975; Dickens & Bell 1976; Sweigart & Mengel 1979;
Norris 1981), in part, because the comparable features of peculiar field stars
were generally consistent with internal processing (Wallerstein & Greenstein,
1964; Wallerstein, 1969). This consensus was eventually challenged by the
observed spread in CNO abundances at early evolutionary stages (Da Costa &
Demarque, 1982; Norris & Freeman, 1979; Hesser & Bell, 1980), evidence of
variations in [Na/Fe] and [Al/Fe] (Cohen, 1978; Peterson, 1980; Cottrell & Da
Costa, 1981; Norris et al., 1981) and the discovery of heavy element (atomic
number greater than $13$) variations in the globular cluster $\omega$ Centauri
(Freeman & Rodgers, 1975).

Closely related to chemical heterogeneity is the so-called second parameter
problem, characterized by the diversity of horizontal branch morphologies
among color-magnitude diagrams (CMDs) of globular clusters with comparable
metallicities (Sandage & Wildey, 1967; van den Bergh, 1967). In some cases,
the horizontal branch is distinctly bimodal (e.g., NGC 2808, Harris 1974),
which led Rood & Crocker (1985) (also see Crocker et al. 1988, as well as
earlier speculation by Cottrell & Da Costa 1981; Norris et al. 1981) to
suggest the presence of multiple populations (MP) of stars with distinct
histories of chemical enrichment. In this scenario, a fraction of members (the
primordial population, using the terminology of Bastian & Lardo 2018) form
similarly to regular halo stars, while the rest (the enriched population) are
influenced by the unique environment of the cluster. The connection between
the horizontal branch morphology and chemical anomalies was first investigated
by Norris (1981); Norris et al. (1981); Smith & Norris (1982a, b, 1983) on the
basis of the observed bimodality in CN absorption among globular clusters with
split horizontal branches. This hypothesis is now firmly established based on
detailed spectroscopic abundance measurements of the horizontal branch members
(Marino et al., 2011).

Over the last three decades, the existence of MP in globular clusters has been
confirmed by detailed spectroscopic analysis of giant (e.g., Dickens et al.
1991; Sneden et al. 1992; Thygesen et al. 2014; Cordero et al. 2014; Marino et
al. 2012; Yong et al. 2008), subgiant and upper main sequence stars (e.g.,
Marino et al. 2016; Gratton et al. 2001; Briley et al. 1996), integrated
spectroscopy (Rennó et al., 2020), and splitting of main sequence (e.g.,
Anderson 1997; Bedin et al. 2004; Piotto et al. 2007; Milone et al. 2013) and
post-main sequence (e.g., Han et al. 2009; Piotto et al. 2012) CMDs. A more
comprehensive summary of available evidence of MP may be found in numerous
reviews of the subject (Bastian & Lardo 2018; Gratton et al. 2004, 2015; Smith
1987; Gratton et al. 2019; Milone & Marino 2022; Dotter 2013 and references
within).

The enriched population is characterized by distinct light element abundances:
most prominently, the overabundance of nitrogen and sodium, and the depletion
of carbon and oxygen, compared to the primordial population. On the other
hand, the abundances of heavy elements generally do not display star-to-star
variations (globular clusters are said to be mono-metallic, in the sense that
$\ce{[Fe/H]}\approx\mathrm{const}$), with the exception of a few anomalous
globular clusters (Marino, 2013), of which $\omega$ Centauri is the most well-
known (Gratton, 1982; Johnson & Pilachowski, 2010; Nitschai et al., 2023). By
contrast, these characteristic abundance patterns are not observed in open
clusters (MacLean et al., 2015; Martell & Smith, 2009; de Silva et al., 2009)
and are rare ($\sim 3\%$ of the $[\ce{Fe/H}]\leq-1$ local sample, Martell et
al. 2011) among field stars (Kraft, 1994; Gratton et al., 2000; Martell &
Grebel, 2010; Carretta et al., 2010), motivating the common assumption that
the handful of known examples were ejected from globular clusters in the past
(Ramírez et al., 2012; Lind et al., 2015).

The physical origin of MP remains largely unexplained (Bastian & Lardo, 2018).
The mono-metallic nature of non-anomalous globular clusters is usually
attributed to the high velocity of supernova ejecta that expel a significant
amount of pristine gas from the cluster, alongside most of the material
enriched with heavy elements, within $\sim 10-35\ \qty{}{Myr}$ of the onset of
star formation (D’Ercole et al., 2008, 2016). Many of the proposed MP theories
fall under one of two broad categories. In multiple generation models (e.g.,
Cottrell & Da Costa 1981; Decressin et al. 2007; Denissenkov & Hartwick 2014),
star formation proceeds in two or more bursts with a sufficient time interval
to allow the earlier generations of stars to build up a reservoir of enriched
gas, from which later generations are assembled. Alternatively, concurrent
formation models (e.g., Bastian et al. 2013; Gieles et al. 2018; Winter &
Clarke 2023) invoke a single generation of coeval stars, in which a fraction
of the members are polluted by stellar ejecta to attain the chemical
abundances of the enriched population. Both approaches suffer from major
shortcomings. Since nuclear processing primarily occurs in massive stars, and
since these stars make a subdominant contribution to the overall mass budget
for commonly assumed initial mass functions (e.g., Kroupa 2001), it is unclear
how a sufficient amount of processed material can be produced in multiple
generation models to assemble the enriched population of stars (enriched stars
make up the majority of members in most globular clusters, Bastian & Lardo
2015). This mass budget problem is easier to resolve in concurrent formation
models, where the enriched population is envisioned to form from pristine gas
prior to the onset of supernovae; however, considerable fine-tuning is needed
to match the required pollution timescales (Salaris & Cassisi, 2014) and to
reproduce the discreetness of populations that arises naturally when multiple
generations are involved.

A potential observational signature of concurrent formation models is the
possible dependence of element abundances on the initial stellar mass
(Ziliotto et al., 2023; Milone, 2023), related to the fact that the magnitude
of chemical enrichment is determined not only by the composition of the
polluted material produced by the primordial generation, but also by the
ability of the enriched generation of stars to accrete it. If the polluted
material is primarily accreted onto circumstellar disks, the accretion
efficiency is expected to be proportional to the surface area of the disk and
its longevity, which, in turn, depends on the mass and evolution of the parent
star (Carpenter et al., 2006; Ribas et al., 2015). If instead the material is
accreted onto the surface of the star, the accretion rate is expected to be
proportional to the squared stellar mass (e.g., in the Bondi accretion
formalism, Bondi 1952). The accreted material may then be diluted by
convective mixing within the star, which also depends on the stellar mass
(Chabrier & Baraffe, 1997). Finally, the density of polluted material may not
be uniform throughout the cluster. For example, in the concurrent formation
model of Bastian et al. (2013), the polluted material is only available in the
core of the cluster. In this case, the accretion efficiency will depend on the
kinematic properties of the stars, which are also related to stellar mass
(Anderson & van der Marel, 2010; Anderson, 2002).

Due to the large distances to globular clusters ($\gtrsim\qty{5}{kpc}$, Harris
1996), high signal-to-noise ratio spectra can only be obtained for the
brightest members that fall within a narrow mass range. Studying the variation
of chemical abundances over a wider mass range inevitably requires the
inference of chemistry from multiband photometry. In general, the photometric
colors of low-mass stars are highly sensitive to the variations in element
abundances due to the dominant molecular chemistry and opacity in low-
temperature atmospheres (Marley et al., 2002; Gerasimov et al., 2022b).
Furthermore, low-mass stars are particularly valuable as chemical tracers,
since molecular absorption is less affected by the non-equilibrium radiation
field (Johnson, 1994), while the interiors of low-mass stars are fully mixed
and undergo minimal nuclear processing (Chabrier & Baraffe, 1997; Gerasimov et
al., 2022a). Sophisticated stellar models are required to take full advantage
of this potential. Accurate simulation of the intricate physics and chemistry
in low-temperature atmospheres remains extremely challenging.

Considerable progress in the photometric analysis of cool stars near the end
of the main sequence ($\gtrsim\qty{0.1}{}$) in globular clusters has been made
since the advent of space-based observations with the Hubble Space Telescope
(HST) and the James Webb Space Telescope (JWST). Recent results include
measurement of the [O/Fe] spread in the globular clusters NGC 6752 (Dotter et
al., 2015), M92 (Ziliotto et al., 2023) and 47 Tucanae (Milone, 2023),
photometric characterization of individual populations in NGC 6752 (Milone et
al., 2019; Gerasimov et al., 2022b), and measurement of $\alpha$-enhancement
in $\omega$ Centauri (Gerasimov et al., 2022a).

The CMDs of globular clusters are expected to extend far beyond the end of the
main sequence into the brown dwarf regime. Unlike stars, brown dwarfs do not
establish energy equilibrium and, instead, undergo long-term cooling, thereby
creating a stellar/substellar gap in the CMD between the faintest main
sequence star and the brightest brown dwarf (Burgasser, 2004; Caiazzo et al.,
2017, 2019). Photometric observations of the substellar sequence in globular
clusters would then not only extend the mass baseline of MP studies, but also
provide an independent constraint on the cluster age from the width of the
gap. The faint magnitudes of globular cluster brown dwarfs pose a major
challenge to their detection. The results of dedicated searches for brown
dwarf candidates in the globular cluster M4 (Dieball et al., 2016, 2019) with
HST remain inconclusive. In our previous work (Gerasimov et al., 2022a), we
estimated the colors and magnitudes of the brown dwarfs in the globular
cluster $\omega$ Centauri and concluded that the substellar sequences of
nearby globular clusters fall within the sensitivity limit of JWST. Since
then, Nardiello et al. (2023) have reported the first tentative discovery of a
brown dwarf in the globular cluster 47 Tucanae, designated BD10. While the
substellar nature of this object remains to be confirmed, we anticipate that
many similar discoveries will follow in the near future from the ongoing Cycle
1 JWST campaigns (Bedin et al., 2021; Caiazzo et al., 2021) and subsequent
cycles.

In this work, we describe a new method to determine the chemical abundances
and other fundamental parameters of globular clusters based on their CMDs from
the subgiant branch to the substellar regime. Our approach involves the
computation of new evolutionary models and model atmospheres with fully self-
consistent chemical abundances. Model isochrones are calculated and fitted to
the observed distribution of photometric colors in an iterative manner. The
computational efficiency of the process is attained by identifying the
components of the stellar models that are most sensitive to particular
elements, and recalculating them only when the abundances of those elements
are updated. We apply our method to the brightest mono-metallic (non-
anomalous) globular cluster 47 Tucanae.

This paper is organized as follows. In Section 2 we describe our compilation
of spectroscopic 47 Tucanae abundances in the literature. The archival HST
photometry that was used in this study is briefly described in Section 3.
Section 4 details the process of calculating a theoretical isochrone for a
given set of chemical abundances, and includes a thorough analysis of
associated systematic errors. Our method of isochrone fitting to the observed
CMD is presented and applied to 47 Tucanae in Section 5. In Section 6, we
derive the mass function of the cluster and predict the anticipated properties
of its substellar members. Our results are discussed and the conclusions are
drawn in Section 7.

##  2 Nominal chemistry

In this study, we adopt spectroscopically inferred chemical abundances of 47
Tucanae as a baseline for comparison against their photometric counterparts,
as well as an initial guess in CMD fitting (Section 5). We refer to this set
of abundances as the nominal composition of the cluster. The nominal
abundances are compiled from three sources: Thygesen et al. (2014), Cordero et
al. (2014) and Marino et al. (2016). The set is similar to that used in Zhou
et al. (2022). Here, we provide a more complete description of its derivation.

Thygesen et al. (2014) list metallicity ([Fe/H]) and the abundances of $26$
other elements, measured for $13$ stars at the tip of the red giant branch
($\qty{3900}{}\lesssim T_{\mathrm{eff}}\lesssim\qty{4300}{}$,
$0.3\lesssim\log_{10}(g)\lesssim 1.5$) in 47 Tucanae. For [Fe/H], [Ti/Fe] and
[Sc/Fe], separate measurements are made for the lines of neutral and ionized
species. Both estimates of [Fe/H] are consistent and treated as independent in
our analysis. Since titanium and scandium have much lower ionization
potentials than iron (Lide, 2004), the neutral lines of these elements may be
affected by unaccounted ultraviolet overionization, particularly prominent at
low effective temperatures and metallicities (Bergemann, 2011; Zhang et al.,
2008; Mashonkina, 2010). As such, we discard the neutral line measurements of
[Ti/Fe] and [Sc/Fe] in Thygesen et al. (2014). We further discard all
measurements without uncertainties.

Cordero et al. (2014) obtained $181$ composition measurements for $164$ unique
red giant and asymptotic giant members of 47 Tucanae ($\log_{10}(g)\lesssim
3$). The dataset includes [Fe/H] and $9$ other elements. Of $164$ observed
stars, $5$ have also been analyzed in Thygesen et al. (2014). Abundance
measurements reported in both references are generally consistent within the
published error bars, with the exception of [Al/Fe], for which the value in
Cordero et al. (2014) exceeds that of Thygesen et al. (2014) by $\sim 2$ and
$\sim 2.5$ sigma for the stars 2618 and 3622 respectively (catalog numbers
from Lee 1977). The difference may be partially caused by the more detailed
modelling of Al lines in Thygesen et al. (2014), following Andrievsky et al.
(2008).

The measurements in Marino et al. (2016) extend our sample to the fainter red
giant and sub-giant members with $3<\log_{10}(g)<4$. The dataset lists
metallicities and abundances of $7$ individual elements for $74$ stars. For
sodium, we adopt the values with non-local thermodynamic equilibrium
corrections based on the curves-of-growth from Lind et al. (2011). Following
the authors’ recommendation, we use uncertainties of $\qty{0.1}{dex}$,
$\qty{0.1}{dex}$ and $\qty{0.15}{dex}$ in [Mg/Fe], [Al/Fe] and [Si/Fe]
respectively, instead of the quoted measurement errors due to the effect of CN
molecular features on the analyzed lines of these elements around
$\qty{0.8}{\micro}$.

The final sample includes measurements from multiple sources for [Fe/H],
[Al/Fe], [Ca/Fe], [Eu/Fe], [La/Fe], [Mg/Fe], [Na/Fe], [Ni/Fe], [O/Fe], [Si/Fe]
and [Ti/Fe]. We carried out a Kolmogorov–Smirnov statistical test for each of
these elements to determine whether the measurements from different sources in
literature are consistent with a shared parent population. The null-hypothesis
(measurements are consistent) was rejected with $99\%$ confidence for [Al/Fe],
[Eu/Fe], [Mg/Fe], [Na/Fe], [Ni/Fe] and [Si/Fe]. Since the spectroscopic
measurements used in this study span a vast range of post-main sequence
evolutionary stages, the discrepancy may reflect a genuine alteration of
surface abundances by the first dredge-up (Salaris et al., 2020), meridional
circulation (Sweigart & Mengel, 1979), thermohaline mixing (Angelou et al.,
2011) or other mechanisms. However, the effect is normally most pronounced in
the products of the CNO cycle (Messenger & Lattanzio, 2002), all of which have
passed the consistency test in our compilation. Therefore, the observed
discrepancies between different literature sources are likely systematic in
nature, exemplified by the aforementioned mismatch in aluminum measurements
between Thygesen et al. (2014) and Cordero et al. (2014).

All spectroscopic composition measurements from the three literature sources
were combined into a single set of chemical abundances, available in Table 3
of Appendix A. We assumed that each abundance measurement ($x_{i}^{(X)}$ for
the $i$-th measurement of element $X$ out of $N^{(X)}$ measurements in total)
is drawn from a normal distribution with the standard deviation composed of
two components added in quadrature: the physical variation in chemistry among
cluster members ($s^{(X)}$, identical for all measurements of $X$) and the
measurement error ($\sigma_{i}^{(X)}$). The physical variation of each element
may then be estimated as in Equation 1:

$\left(s^{(X)}\right)^{2}=\frac{N^{(X)}\mathrm{Var}\left(x_{i}^{(X)}\right)}{N^{(X)}-1.5}-\left\langle\left(\sigma_{i}^{(X)}\right)^{2}\right\rangle$
(1)

Since we are ultimately interested in the standard deviation of the physical
spread, we evaluated the unbiased sample variances assuming $N^{(X)}-1.5$
degrees of freedom (Brugger, 1969). In cases where the average measurement
error exceeds the observed scatter in the data, $s^{(X)}$ cannot be estimated,
and the data may be considered consistent with lack of star-to-star variations
in $X$. The mean abundance of each element and its uncertainty were calculated
using the square reciprocal measurement errors as statistical weights with
$s^{(X)}$ added to the measurement error in quadrature if available. Finally,
we note that the estimates of $s^{(X)}$ for the six elements that did not pass
the Kolmogorov–Smirnov test are likely biased by systematic errors that may
have been unaccounted for in the quoted measurement uncertainties.

Since the helium abundance cannot be measured spectroscopically at low
effective temperatures, the nominal helium mass fraction of 47 Tucanae was
adopted as $Y=0.25$, following the isochrone fit of eclipsing binary members
in Thompson et al. (2020).

##  3 Archival photometry

Throughout this paper, all mentions of the optical bands F606W and F814W will
implicitly refer to the wide field channel of the Advanced Camera for Surveys
(ACS) on HST. Likewise, all mentions of the near infrared bands F110W and
F160W will refer to the near infrared channel of the Wide Field Camera 3
(WFC3) on HST. Finally, the infrared bands F150W2 and F322W2 will refer to the
Near Infrared Camera (NIRCam) on JWST.

Our analysis is based on archival HST photometry (GO-11677, PI: H. Richer) of
47 Tucanae presented in Kalirai et al. (2012). The primary field was imaged in
the F606W and F814W bands. The field spans $\sim\qty{5}{arcmin}$ and is
located at $\qty{6.7}{arcmin}$ ($2.1$ half-light radii, Harris 1996) from the
center of the cluster. Parallel fields were imaged in the F110W and F160W
bands, as well as additional bands of the ultraviolet and visible light
channel that are not used in this study. The parallel fields span a $250$
degree arc, centered at the primary field with the radius of
$\sim\qty{6}{arcmin}$ and facing away from the center of the cluster (see
Figure 2 of Kalirai et al. 2012). This work is based on the observations of
the primary field and one of the parallel fields with the largest cumulative
exposure time (field 13, referred to as the “stare”).

The archival photometry is contaminated by the members of the Small Magellanic
Cloud (SMC), centered $\sim 2^{\circ}$ southeast of 47 Tucanae. Since the SMC
is more distant than 47 Tucanae by over an order of magnitude, the main
sequences of the two objects are well-separated in the optical CMD and have a
small overlap in the near infrared CMD (see Figure 6). The overlap region is
excluded from our analysis as detailed in Section 6. A more thorough reduction
of the archival data with proper motion cleaning is deferred to a future
study.

All magnitudes quoted in this paper are VEGAMAG.

##  4 Model isochrones

We derive the properties of 47 Tucanae by comparing theoretical model
isochrones to the observed CMD of the cluster. Due to the presence of MP, an
accurate model of the CMD requires multiple theoretical isochrones that
capture the observed spread in photometric colors. We therefore aim to find
three isochrones with distinct chemical compositions that approximately trace
out the blue tail, the red tail and the ridgeline of the CMD. We adopt the
convention of identifying the blue and red tails of the distribution based on
the F110W – F160W near infrared color of the lower main sequence (the optical
colors are, in fact, inverted compared to their near infrared counterparts).

We search for the desired isochrones through iterative perturbations of the
nominal chemical composition. Both the final isochrones and the intermediate
isochrones at each iteration are calculated using an improved version of the
method, originally developed for our previous study of the globular cluster
$\omega$ Centauri (Gerasimov et al., 2022a). In this section, we describe the
general process of calculating a theoretical isochrone for a given chemical
composition. Our approach to isochrone fitting and the treatment of multiple
populations in the cluster are presented in Section 5.

###  4.1 Evolutionary models

Stellar evolution is simulated using the MESA (Modules for Experiments in
Stellar Astrophysics) code (Paxton et al., 2011, 2013, 2015, 2018, 2019),
version 15140. The evolution begins at the pre-main sequence (PMS) phase and
is allowed to proceed either until the age of $\qty{13.5}{Gyr}$ or until the
tip of the subgiant branch, whichever occurs sooner. For our purposes, we
define the tip of the subgiant branch as the point, at which the innermost
shell that satisfies $\varepsilon_{\mathrm{nuc}}>\qty{e3}{\per\per}$ begins to
fall outside the central $10\%$ of the stellar mass
($\varepsilon_{\mathrm{nuc}}$ is the specific nuclear energy release rate).
This criterion was empirically determined to be a reliable indicator of the
hydrogen shell ignition that characterizes the onset of the red giant branch
(Gamow & Keller, 1945). Evolving the model further into the red giant branch
incurs a larger computational cost due to the rapid variation of surface
parameters with age, and is beyond the scope of this study; however, see
Salaris et al. (2020).

The initial stellar masses are sampled on an adaptive grid that guarantees the
difference in luminosities and effective temperatures between adjacent grid
masses of $\left|\Delta\log_{10}{L}\right|<\qty{0.12}{dex}$ and $\left|\Delta
T_{\mathrm{eff}}\right|<\qty{120}{}$, respectively, at the checkpoint ages
between $\qty{10}{Gyr}$ and $\qty{13.5}{Gyr}$ in steps of $\qty{0.5}{Gyr}$
(the typical average differences at $\qty{13.5}{Gyr}$ are
$\langle\left|\Delta\log_{10}{L}\right|\rangle\approx\qty{0.07}{dex}$ and
$\langle\left|\Delta T_{\mathrm{eff}}\right|\rangle\approx\qty{70}{}$). The
lowest mass in the grid ($\sim\qty{0.04}{}-\qty{0.05}{}$) is chosen to attain
an evolved $T_{\mathrm{eff}}{}$ of $\lesssim\qty{700}{}$ at $\qty{13.5}{Gyr}$,
while the highest mass ($\sim\qty{0.9}{}$) is set by requiring that the star
takes at least $\qty{10}{Gyr}$ to reach the tip of the subgiant branch.

Figure 1: Effective temperatures and surface gravities (Kiel diagram) of the
model atmospheres, calculated in this study for the final ridgeline isochrone
of 47 Tucanae. The PHOENIX models are color-coded by their convergence
patterns, as detailed in Section 4.4. The evolutionary tracks of selected MESA
models in the $T_{\mathrm{eff}}{}$ – $\log_{10}(g){}$ space are shown for
comparison, as well as the final isochrone at $\qty{11.5}{Gyr}$. The tracks
are labeled by their initial stellar masses in $\qty{}{}$.

PMS stars at zero age are initialized with uniform chemical abundances and
central temperature of $T_{c}=\qty{5e5}{}$. This choice follows Choi et al.
(2016) and Gerasimov et al. (2022a), and ensures that the stars of all
considered initial masses are found within the boundary condition tables
(Section 4.2) at the completion of the main PMS relaxation routine in MESA.
The additional PMS relaxation for $>\qty{0.3}{}$ stars until the formation of
the radiative core (pre_ms_relax_to_start_radiative_core) was disabled in all
of our evolutionary models to provide a consistent definition of zero age for
all calculated evolutionary tracks. This feature was introduced in recent
versions of the code to improve the robustness of massive star
calculations111https://github.com/MESAHub/mesa/issues/340, and does not
significantly affect the final evolutionary states of the objects considered
in this study.

The interior convective mixing length (in the formalism of mixing length
theory, MLT; Böhm-Vitense 1958) in all models was set to the solar-calibrated
value of $\alpha_{\mathrm{MLT}}=1.82$ scale heights (Salaris & Cassisi, 2015;
Choi et al., 2016). The effect of other choices for this parameter is
discussed in Section 5. The input and output files for the evolutionary models
produced in this study are made available on our
website222http://romanger.com/models.html and
https://doi.org/10.5281/zenodo.10016008 (catalog Zenodo).

Figure 2: Synthetic spectra of selected model atmospheres at
$\log_{10}(g){}=5$, calculated in this study for the ridgeline isochrone of 47
Tucanae. The important absorption bands of $\ce{CH_{4}}$, $\ce{CO}$,
$\ce{CO_{2}}$, $\ce{H_{2}O}$ and $\ce{NH_{3}}$ are labeled with ragged black
lines that represent the variation of absorption strength within the band for
the $T_{\mathrm{eff}}{}=\qty{800}{}$ model. The central wavelengths of
prominent $\ce{TiO}$ bands, Na and K atomic lines and the peak of the CH_4
$\nu_{3}$ (asymmetric stretching vibration band) Q-branch are shown with
vertical dotted lines. The approximate wavelength range of strong collision-
induced H_2 absorption (CIA) is shown with the horizontal dashed line. The
spectra have been convolved with a $\qty{10}{}$ Gaussian kernel for clarity.

###  4.2 Boundary conditions

At each evolutionary step, energy conservation requires the pressure-
temperature profile to be consistent with the luminosity output of the star.
Evaluating this condition near the surface of low-mass stars is challenging
due to the complex relationship between the pressure-temperature profile and
the outgoing energy flux, caused by low-temperature phenomena such as non-grey
molecular opacity, cloud formation and convective instabilities. This issue is
particularly prominent at $T_{\mathrm{eff}}{}\lesssim\qty{4500}{}$, as the
structure of the stellar atmosphere begins to significantly deviate from the
Eddington approximation (Mihalas, 1978; Burrows et al., 1989; Dorman et al.,
1989; Zhou et al., 2022). More accurate atmospheric pressure-temperature
profiles may be extracted from model atmospheres, precomputed for the entire
range of surface parameters that may be encountered by the star during its
evolution.

Following Chabrier & Baraffe (1997); Choi et al. (2016); Phillips et al.
(2020) and our earlier study (Gerasimov et al., 2022a), we calculated a grid
of model atmospheres (Section 4.3) for each theoretical isochrone, spanning
$2\leq\log_{10}(g){}\leq 6$ for $\qty{2500}{}\leq
T_{\mathrm{eff}}{}\leq\qty{7500}{}$, and $4\leq\log_{10}(g){}\leq 6$ for
$\qty{500}{}\leq T_{\mathrm{eff}}{}\leq\qty{2400}{}$, with steps of
$\qty{100}{}$ in $T_{\mathrm{eff}}{}$ and $0.5$ in $\log_{10}(g){}$. Models
with $\log_{10}(g){}<4\ \&\ T_{\mathrm{eff}}{}<\qty{2500}{}$ were not
required, since such low gravities are only encountered during the early
evolution ($<\qty{5}{Myr}$) and the subgiant phase, both of which are
characterized by higher effective temperatures
($T_{\mathrm{eff}}{}>\qty{2500}{}$). For some of the calculated isochrones, a
small number of model atmospheres could not be computed due to convergence
issues (Section 4.4). The corresponding empty grid points were then filled in
using linear Delaunay triangulation (Barber et al., 2013). The temperatures
and gravities of the atmosphere models calculated for the final ridgeline
isochrone are plotted in Figure 1 alongside the $T_{\mathrm{eff}}{}$ –
$\log_{10}(g){}$ tracks of selected evolutionary models.

The model atmospheres were converted to tables of gas pressure and temperature
at a prescribed depth and provided to MESA as the outer boundary conditions at
each evolutionary step and for each iteration in the solution of the stellar
structure equations. In general, specifying the boundary conditions at larger
Rosseland optical depths ($\tau$) is preferred, as it minimizes the
discontinuity in opacity between the atmosphere and interior, and ensures
adiabatic behavior at the boundary (Chabrier & Baraffe, 1997, 2000; Spada et
al., 2017). However, the reduced atmospheric opacity in more massive stars
shifts the boundary to larger physical depths, where the non-ideal gas
behavior may be significant (Burke et al., 2004), and is not accounted for in
our high-$T_{\mathrm{eff}}{}$ model atmospheres. We therefore employed two
distinct atmosphere-interior coupling regimes with boundary conditions at
$\tau=100$ for $M<\qty{0.5}{}$, and at $T(\tau)=T_{\mathrm{eff}}{}$
(photosphere) for $M>\qty{0.6}{}$, where $T(\tau)$ is the temperature profile
of the star and $M$ is the initial stellar mass. At intermediate stellar
masses, $\qty{0.5}{}\leq M\leq\qty{0.6}{}$, a linear ramp between the two
regimes was applied. The range of transition was chosen to approximately
coincide with the maximum value of $dT(\tau=100)/dM$ for main sequence stars.
The implications of this choice on the synthetic photometry are discussed in
Section 4.6.

When the range of pre-tabulated boundary conditions is exceeded, MESA
implements a fail-safe and falls back on the Eddington approximation. Near the
edges of the table, a smooth blending between the on-table and off-table
boundary conditions is carried out, which may result in numerical
instabilities at low $T_{\mathrm{eff}}{}$ due to the extreme inaccuracy of the
Eddington formula. Since our models, by design, do not leave the table range
during regular evolution, we modified the source code of MESA to disable the
off-table blending once the PMS phase has been completed. We also updated the
code to interpolate the tables linearly rather than bicubically to address
unphysical temperatures and pressures resulting from spline overshoots (Fried
& Zietz, 1973) in the vicinity of poorly converged model atmospheres (Section
4.4). The calculated boundary condition tables and the patch file for the MESA
codebase are made available on our website333http://romanger.com/models.html
and https://doi.org/10.5281/zenodo.10016008 (catalog Zenodo).

###  4.3 Model atmospheres

Model atmospheres were calculated using the PHOENIX code (Hauschildt et al.,
1997), version 15.05.05D (Allard et al., 2011, 2012; Gerasimov et al., 2020)
at $T_{\mathrm{eff}}{}\leq\qty{5000}{}$; and the ATLAS code (Kurucz, 1970),
version 9 (Kurucz, 2005, 2014) at $T_{\mathrm{eff}}{}>\qty{5000}{}$. The use
of ATLAS for “warm” atmospheres is advantageous due to its superior
computational efficiency (see Section 4.5), attained by virtue of opacity
sampling from pre-tabulated opacity distribution functions (ODFs; Kurucz et
al. 1974; Carbon 1984), simplified chemical equilibrium (only gaseous species
are considered, molecule-molecule interactions are ignored) and the ideal
equation of state. The ODFs for the chemical composition of each isochrone
were computed for $\qty{9}{\nano}\lesssim\lambda\lesssim\qty{160}{\micro}$ at
$57$ temperatures between $\approx\qty{2e3}{}$ and $\approx\qty{2e5}{}$ using
the satellite package DFSYNHTE (Castelli, 2005; Castelli & Kurucz, 2003).
ATLAS atmospheres are stratified into $72$ plane-parallel layers, spanning
$\sim${10}^{-7}$<\tau<$100$$ with logarithmic spacing. The synthetic spectra
for each model atmosphere were calculated using the SYNTHE code (Kurucz &
Avrett, 1981) within $\qty{0.1}{\micro}<\lambda<\qty{5}{\micro}$, at the
constant resolution444SYNTHE treats all line profiles as symmetric with
respect to the nearest point on the wavelength grid. For this reason, high
resolutions of order $\lambda/\Delta\lambda\sim 10^{5}-10^{6}$ are typically
recommended even when not required for intended science (Kurucz, 2014;
Sbordone & Bonifacio, 2005). of
$\lambda/\Delta\lambda\sim$6\text{\times}{10}^{5}$$. All three codes are
packaged in the open source Python dispatcher
BasicATLAS555https://github.com/Roman-UCSD/BasicATLAS (Larkin et al., 2023),
alongside a suite of internal consistency checks and a synthetic photometry
calculator.

PHOENIX allows for a more accurate treatment of low-$T_{\mathrm{eff}}{}$
atmospheres through direct sampling of opacity at the wavelengths of interest,
a more complete chemical reaction network (including condensation, see Allard
et al. 2001 for a review), a comprehensive molecular line database (most
importantly, H_2O lines from Barber et al. 2006 and CH_4 lines from Brown 2005
and Wenger & Champion 1998; Homeier et al. 2003; other included lines are
listed in Table 3 of Gerasimov et al. 2022a), and a non-ideal equation of
state (Saumon et al., 1995), among other features. The inclusion of condensate
species (clouds) in the chemical equilibrium and the opacity model is
necessary to reproduce the observed reddening of photometric colors at
$T_{\mathrm{eff}}\lesssim\qty{3000}{}$ compared to gas-only models (Allard et
al., 2001; Lunine et al., 1986). Our setup of PHOENIX also includes the Allard
& Homeier treatment of gravitational settling for condensate species (Allard
et al., 2012; Helling et al., 2008) that gradually removes clouds from the
spectrum-forming region of the atmosphere at $T_{\mathrm{eff}}\ll\qty{2000}{}$
(Tsuji & Nakajima, 2003). For computational efficiency, we only consider
gravitational settling at $T_{\mathrm{eff}}<\qty{2500}{}$ and, otherwise,
assume that the condensate species remain at chemical equilibrium.

Our PHOENIX models are stratified into $250$ layers when gravitational
settling is used and $128$ layers otherwise. Unlike ATLAS, the atmospheric
layers in PHOENIX are spherical and defined by the optical depth at a fixed
wavelength ($\tau_{\lambda}$, where $\lambda=\qty{1.275}{\micro}$) instead of
the Rosseland mean666Gerasimov et al. (2022a) have overlooked this distinction
in definitions and, therefore, the boundary condition tables used in that
study may be less accurate, especially at very low $T_{\mathrm{eff}}{}$, where
the opacity at $\lambda=\qty{1.275}{\micro}$ is no longer representative of
the average across all wavelengths.. The bottom of the atmosphere is set to
$\tau_{\lambda}=${10}^{3}$$ for the models without gravitational settling, but
reduced to $\tau_{\lambda}=316$ when settling is enabled to avoid the
numerical issues associated with settling calculations at large depths. The
top of the atmosphere is defined by the gas pressure, rather than optical
depth, set to $\qty{e-3}{dyn\per^{2}}$. The optical depth grid is
approximately logarithmic. For all models considered in this study, the
atmospheres at the final evolutionary state are negligibly thin compared to
the stellar radius (the $\tau=100$ layer is always found within the outer
$0.1\%$ of the photospheric radius, calculated at
$T(\tau)=T_{\mathrm{eff}}{}$). Therefore, the effect of spherical geometry is
expected to be subdominant, allowing the use of approximate atmospheric radii
from literature (Baraffe et al., 2015) instead of introducing additional
dimensions to our model atmosphere grids. The wavelength sampling and spectral
synthesis for all PHOENIX models are carried out within
$\qty{1}{}<\lambda<\qty{1}{\milli}$ with a median resolution of
$\lambda/\Delta\lambda\sim$18\,250$$ between $\qty{0.5}{\micro}$ and
$\qty{3}{\micro}$. A few examples of the synthetic spectra for the chemical
abundances of the ridgeline isochrone and $\log_{10}(g){}=5.0$ are provided in
Figure 2. The atmosphere models calculated in this study are made available
online777http://romanger.com/models.html.

We assume that the chemical composition of the atmosphere is unaffected by
stellar evolution and is identical to the initial composition of the PMS star,
with the exception of [Li/Fe] that is reduced by $\qty{3}{dex}$ in all PHOENIX
models compared to the PMS abundance. This assumption is justified by the fact
that low-mass stars ($\lesssim\qty{0.6}{}$) undergo minimal nuclear processing
(with the exception of lithium burning), while their higher-mass counterparts
establish radiative zones at young ages, thereby preventing the nuclear
burning products from reaching the atmosphere. The behavior of [Li/Fe] in the
interior and the atmosphere of the star as a function of mass is explored in
detail in Gerasimov et al. (2022a).

###  4.4 Model convergence

In both model atmospheres and the time steps of evolutionary models, the
equations of input physics are solved using iterative methods. At the high
effective temperatures of ATLAS atmospheres, the parameter space maintains
approximate local linearity (Gustafsson et al., 2008), allowing for fast and
reliable convergence onto a self-consistent solution. We carry out ATLAS
iterations in batches of $15$ until the maximum flux error and the maximum
flux derivative error across all layers drop below $1\%$ and $10\%$,
respectively (following Sbordone & Bonifacio 2005; Mészáros et al. 2012; see
the appendix in Larkin et al. 2023 for details), or until no further progress
can be made. These convergence standards have been met by nearly all ATLAS
atmospheres computed in this study, with the exception of a few
$T_{\mathrm{eff}}{}>\qty{6500}{}$ models that generally fall outside the
region of parameter space, explored by the calculated evolutionary tracks
(Figure 1).

The evolution of convergence parameters across iterations in
low-$T_{\mathrm{eff}}{}$ PHOENIX atmospheres is more involved. In models with
subdominant cloud settling ($T_{\mathrm{eff}}{}\gtrsim\qty{2000}{}$), the flux
errors generally decrease with every iteration until the model arrives at a
stable solution with the typical maximum and average flux errors across all
radiative layers of $\sim 3\%$ and $0.5\%$, correspondingly. In this work, we
refer to this convergence pattern as stable, and the models that exhibit this
pattern are shown in green in Figure 1 for the ridgeline isochrone. A small
number of models with stable convergence may occasionally encounter iterations
with ill-conditioned temperature corrections that break the trend in
convergence and increase the flux errors by over an order of magnitude,
effectively restarting the computation. In those cases, the iteration with the
minimum average error, from which we derive the final synthetic spectrum, is
still expected to produce a reliable result; however, the gain in accuracy
from an increased number of iterations may be limited.

At lower $T_{\mathrm{eff}}{}$, the effects of condensate settling become more
prominent, resulting in a far more complex relationship between the model
parameters that restricts the effectiveness of the temperature correction
scheme. In this regime, the flux errors typically oscillate between high and
low values across iterations. While the flux errors of the best iteration are
usually comparable to those of the atmospheres with stable convergence, the
final solutions to the structure equations may lack uniqueness (i.e. for some
$T_{\mathrm{eff}}{}$ and $\log_{10}(g){}$ there may be multiple atmosphere
structures with similar flux errors but vastly different emergent spectra).
Furthermore, the pronounced sensitivity of the model to temperature
corrections casts doubt on the usefulness of the radiative flux errors as a
diagnostic of self-consistency. We refer to the convergence pattern in this
temperature regime as irregular. Some of the PHOENIX atmospheres with
irregular convergence had to be excluded from the grid due to the rapid growth
of the oscillation amplitude in the temperature structure between iterations.
For the rest of models with this convergence pattern, the solution with the
lowest average radiative flux error was added to the grid at the points
identified in Figure 1 with black markers for the ridgeline isochrone.

A small number of models displayed a convergence pattern, intermediate between
that of stable convergence and irregular convergence, which we refer to as
limited convergence. In those cases, the model may exhibit stable convergence
for the first few iterations, but then transition into the irregular
convergence pattern before a solution with satisfactory flux errors can be
reached. In some cases, the transition between the two convergence patterns
may occur multiple times over the course of $\sim 100$ iterations. The models
with limited convergence are shown in Figure 1 with red markers for the
ridgeline isochrone. At very low temperatures
($T_{\mathrm{eff}}{}\lesssim\qty{1200}{}$) and high gravities
($\log_{10}(g){}\gtrsim 5.5$), the condensate species are almost completely
removed from the spectrum-forming region of the atmosphere, restoring the
stable convergence pattern (see Figure 1).

For most PHOENIX models computed in this study, we used the nearest atmosphere
structure from the NEXTGEN model grid (Hauschildt et al., 1999; Allard et al.,
2014) as the initial guess in the solver, with the exception of models with
particularly poor convergence that were recalculated, using the nearest well-
converged atmospheres from our grid as the initial models instead. In
principle, the accuracy of our low-temperature atmospheres may be improved by
calculating the models in small batches with progressively decreasing
$T_{\mathrm{eff}}{}$ and $\log_{10}(g){}$, and using the atmospheres from the
preceding batch as initial guesses. In fact, we adopted this approach for the
ATLAS models. However, the high computational demand of PHOENIX (Section 4.5)
makes the un-parallelized computation extremely time-consuming. Furthermore,
it is unclear how the potential gain in numerical precision would compare to
the systematic errors in the input physics and grid interpolation errors. We
therefore chose to focus on identifying the outlier models based on their
synthetic photometry and excluding them from the isochrone calculation
instead, as described in Section 4.7.

Nearly all evolutionary MESA models in our grids have converged at every time
step within the “gold tolerances” (Paxton et al., 2019), with the exception of
a few lowest-mass models ($\lesssim\qty{0.05}{}$), where the tolerances were
relaxed to their standard values to allow evolution over the discontinuities
in $dT(\tau=100)/d\log_{10}(g){}$, caused by the sparse gravity sampling in
the boundary condition tables.

###  4.5 Computational demand of model atmospheres

All model atmospheres used in this work were calculated on the Bridges-2
supercomputer (Brown et al., 2021), operated by the Pittsburgh Supercomputing
Center. Our PHOENIX models were calculated to $\sim 100$ iterations, requiring
between $\sim 0.15$ processor hours per iteration at the highest temperatures
to $\sim 1.5$ processor hours at the lowest. The high memory demand of PHOENIX
($\gtrsim\qty{4}{GB}$ per model) necessitated requesting at least $2$ (and
occasionally $3$) service units for each processor hour.

The computational demand of ATLAS models is dominated by the spectral
synthesis, carried out by the SYNTHE package. The calculation of the emergent
spectra for most ATLAS models required between $1.5$ and $3.5$ processor hours
per spectrum (and the same number of service units). For comparison, the
structure calculations only took $\approx 0.005$ processor hours per model for
all iterations.

The computational demand of a complete atmosphere grid for each chemical
composition is approximately $5\text{\times}{10}^{4}$ service units, with
nearly $98\%$ taken by PHOENIX calculations. This estimate does not include
recalculation of failed models, ODF calculations, calculations of partial
pressure tables for each chemical composition used in PHOENIX, and
evolutionary calculations with MESA. This estimate only applies to the final
isochrones (ridgeline, blue tail and red tail) presented in Section 5: to
calculate the intermediate isochrones used during the fitting process, we took
advantage of various optimizations described in Section 5.

###  4.6 Isochrone stitching

The process of calculating theoretical isochrones described earlier in this
section required a choice of three “stitching points”, where distinct
modelling setups were blended together: the effective temperature where ATLAS
models are replaced with PHOENIX models ($\qty{5000}{}$), the effective
temperature where the gravitational settling in PHOENIX models was disabled
($\qty{2500}{}$) and the initial stellar mass, where the $\tau=100$ boundary
condition tables were replaced with $T(\tau)=T_{\mathrm{eff}}{}$ (smooth ramp
between $\qty{0.5}{}$ and $\qty{0.6}{}$). Here, we briefly review the
implications of those choices for synthetic photometry.

Our choice of the ATLAS/PHOENIX transition temperature ($\qty{5000}{}$) is
higher than that used in our previous work ($\qty{4000}{}$, Gerasimov et al.
2022a). For comparison, the published grid of ATLAS models with the updated
molecular opacities (Castelli & Kurucz, 2003) reaches $\qty{3500}{}$, although
the low-temperature extension has known issues (Plez, 2011). We calculated an
additional test set of ATLAS models for the parameters of the ridgeline
isochrone at $\qty{4000}{}\leq T_{\mathrm{eff}}{}\leq\qty{5000}{}$, as well as
an additional grid of evolutionary models with the updated boundary
conditions. We found the synthetic photometry at $\qty{11.5}{Gyr}$ of the
regular ridgeline isochrone and its altered version to be indistinguishable at
$T_{\mathrm{eff}}=\qty{5000}{}$, but rapidly diverging at lower temperatures
to attain the difference of $\approx\qty{0.025}{mag}$ in F606W – F814W and
F150W2 – F322W2 at $T_{\mathrm{eff}}{}=\qty{4000}{}$. The difference
originates almost entirely from the synthetic spectra with only a minor
contribution from the boundary conditions.

Our choice of the minimum $T_{\mathrm{eff}}{}$ without gravitational settling
($\qty{2500}{}$) is, on the other hand, lower than that employed in our
previous study of $\omega$ Centauri ($\qty{3000}{}$, Gerasimov et al. 2022a),
which allowed us to significantly reduce the computational demand. We
calculated a test grid of PHOENIX atmospheres for the ridgeline isochrone with
$\qty{2500}{}\leq T_{\mathrm{eff}}{}\leq\qty{3000}{}$ and enabled settling of
condensates, as well as a test set of corresponding evolutionary models. We
found that at $T_{\mathrm{eff}}{}=\qty{2500}{}$ and at the age of
$\qty{11.5}{Gyr}$, the effect of gravitational settling on the same colors as
before is $\approx\qty{0.01}{mag}$. In the optical regime, both the boundary
conditions and the synthetic spectra make comparable contributions to the
observed difference, while the latter dominate in the infrared. At
$T_{\mathrm{eff}}{}>\qty{2500}{}$, the effect of gravitational settling
decreases rapidly in the infrared, but remains approximately constant in the
optical up to the warm end of the test grid at $\qty{3000}{}$, primarily due
to the effect of the updated boundary conditions on stellar evolution.

Lastly, we examine the effect of the transition between the atmosphere-
interior coupling schemes by disabling the linear ramp at $\qty{0.5}{}\leq
M\leq\qty{0.6}{}$. The effect on synthetic photometry at $\qty{11.5}{Gyr}$ is
most prominent in F606W – F814W when the range of $T(\tau)=T_{\mathrm{eff}}{}$
boundary condition tables is extended down to $\qty{0.5}{}$, resulting in the
discontinuity of $\lesssim\qty{0.01}{mag}$ between the two regimes. The
discontinuity is less prominent when the $\tau=100$ tables are extended to
$\qty{0.6}{}$ or when other photometric colors are considered.

In summary, our choice of “stitching points” in the isochrone is not expected
to introduce errors larger than $\qty{0.01}{mag}$ in any of the photometric
colors considered in this study; however, it appears that the optimal choice
of the transition between the two atmosphere-interior coupling schemes falls
at somewhat larger masses, and a small gain in accuracy may be attained by
allowing gravitational settling of condensates at
$T_{\mathrm{eff}}>\qty{2500}{}$.

Figure 3: Comparison of synthetic photometry at 11.5Gyr and
$T_{\mathrm{eff}}{}\leq\qty{2350}{}$ for the blue tail isochrone, calculated
using the two methods described in Section 4.7: atmosphere hammering (black)
and linear interpolation of the bolometric corrections (red). The outlier
models were identified by visual inspection in each isochrone and excluded
from the “masked” result.

###  4.7 Synthetic photometry and hammering

The preliminary synthetic photometry for the calculated isochrones was
obtained by applying the reddening law to the synthetic spectra, evaluating
the bolometric corrections for each model atmosphere in the bandpasses of
interest, and interpolating the result to $T_{\mathrm{eff}}{}$ and
$\log_{10}(g){}$, predicted by the evolutionary model at the required initial
mass and age. The resulting magnitudes were corrected by the distance modulus
of $\qty{13.2418\pm 0.0625}{mag}$ (Chen et al., 2018). To calculate the
bolometric corrections, we used the synphot() routine of the BasicATLAS
package (Larkin et al., 2023). The process uses the reddening law from
Fitzpatrick & Massa (2007), parameterized only by the optical reddening,
$E(B-V)$, and assuming the total-to-selective extinction ratio of $R_{V}=3.1$.

The strong dependency of low-$T_{\mathrm{eff}}{}$ spectra on surface
parameters and the highly non-linear behavior of cool atmospheres pose a
challenge to the interpolation of bolometric corrections. Linear interpolation
of a sparsely sampled grid leads to large unphysical discontinuities in the
derivatives, and is known to introduce noticeable errors even at relatively
high temperatures (Mészáros & Allende Prieto, 2013; Gustafsson et al., 2008).
On the other hand, higher-order interpolation methods are less robust against
outliers (Fried & Zietz, 1973), such as those introduced by the models with
poor convergence. Specialized interpolation methods for atmosphere grids with
non-convergent models (e.g. Nendwich et al. 2004) generally require a means to
reliably identify the outliers, which is not straightforward in the case of
irregular convergence, as discussed in Section 4.4.

To address these issues, we decided to avoid interpolating the grid at low
temperatures altogether and, instead, to calculate additional PHOENIX
atmospheres at $\qty{600}{}\leq T_{\mathrm{eff}}{}<\qty{2400}{}$ in
$\qty{50}{}$ steps, with $\log_{10}(g){}$ set to the values predicted by the
isochrone for each effective temperature at the target age. We refer to this
process as atmosphere hammering. In addition to removing the need to
interpolate bolometric corrections at low temperatures, hammering serves two
other purposes. First, it reduces the effective number of dimensions in the
atmosphere grid from $2$ to $1$ (since the hammering models obey a known
$T_{\mathrm{eff}}{}$ – $\log_{10}(g){}$ relationship), allowing the outlier
models to be easily identified, e.g., by their placement in the color-
magnitude space. Second, it allows us to derive the synthetic photometry at
$T_{\mathrm{eff}}{}<\qty{2500}{}$ from higher-gravity model atmospheres that
have better convergence (e.g., from Figure 1, to calculate the bolometric
corrections at the $T_{\mathrm{eff}}{}=\qty{2000}{}$ point on the
$\qty{11.5}{Gyr}$ isochrone, one would require both $\log_{10}(g){}=5.0$ and
$\log_{10}(g){}=5.5$ model atmospheres for the interpolation method, but only
one $\log_{10}(g){}\approx 5.5$ model for the hammering method). The major
drawbacks of the hammering method are the added computational cost and the
resulting commitment to a particular target age, since the hammering models
fix the $T_{\mathrm{eff}}{}$ – $\log_{10}(g){}$ relationship.

The linear interpolation and hammering methods of calculating synthetic
photometry at $T_{\mathrm{eff}}{}\leq\qty{2350}{}$ are compared in Figure 3.
The figure demonstrates that at $T_{\mathrm{eff}}{}\lesssim\qty{1000}{}$ (the
brown dwarf regime), linear interpolation may be off by as much as
$\qty{2}{mag}$ in the optical regime and $\gtrsim\qty{0.5}{mag}$ in the near
infrared. It is also apparent that the two methods converge at
$T_{\mathrm{eff}}{}\gtrsim\qty{2000}{}$, suggesting that linear interpolation
of bolometric corrections is valid across the main sequence. In this study, we
use linear interpolation for all models with
$T_{\mathrm{eff}}{}>\qty{2350}{}$.

We note that while the synthetic photometry inferred from atmosphere hammering
does not rely on interpolated bolometric corrections, it is still based on the
surface parameters predicted by the evolutionary track, which, in turn, was
calculated using interpolated boundary conditions. The effect of boundary
conditions on synthetic photometry is subdominant compared to that of
bolometric corrections; however, some interpolation errors may remain at low
$T_{\mathrm{eff}}{}$. To avoid interpolation of model atmospheres entirely,
one must calculate new atmospheres for every iteration of the evolutionary
code at every time step, which would increase the total computational demand
by over an order of magnitude.

##  5 Isochrone fitting

Figure 4: Top: ridgeline isochrones of 47 Tucanae at different ages,
overplotted on the distribution of observed colors in optical (left) and near
infrared (right). Middle: comparison of the ridgeline isochrone at
$\qty{11.5}{Gyr}$ with the nominal [Fe/H] and $Y$ (red) to the helium-rich
(blue) and metal-rich (black) isochrones at the same age. Bottom: effect of
the mixing length parameter on the ridgeline isochrone at the age of
$\qty{11.5}{Gyr}$. The cyan and yellow regions of the color-magnitude diagram
subtend the $1\sigma$ and $2\sigma$ contours in the data distribution,
respectively. $E(B-V)=\qty{0.04}{mag}$ for all shown isochrones.

###  5.1 Main sequence turn-off and subgiant branch

The high effective temperatures of the main sequence turn-off and the subgiant
branch in 47 Tucanae ($T_{\mathrm{eff}}{}\gtrsim\qty{5000}{}$) suppress
molecular formation in the atmosphere and minimize the impact of chemical
abundances on the observed photometry. The notable exceptions are the overall
metallicity, [Fe/H], and the helium mass fraction, $Y$, as they have a
significant effect on the mean molecular weight and opacity of the interior
that, in turn, guide the evolutionary track of the star (Salaris & Cassisi,
2005, Ch. 5.5). These high-mass stars are also sensitive to the adopted energy
transport parameters (in particular, the mixing length,
$\alpha_{\mathrm{MLT}}$) due to the emergence of an outer convective zone
(Joyce & Chaboyer, 2018). To reduce the degree of degeneracy in the
determination of chemical abundances, it is advantageous to constrain as many
parameters as possible from the main sequence turn-off and the subgiant branch
before analyzing the lower main sequence. We used this part of the CMD to
assess our choice of $\alpha_{\mathrm{MLT}}=1.82$, fix the cluster age and the
optical reddening, $E(B-V)$, as well as to evaluate the accuracy of the helium
mass fraction and metallicity in the nominal chemical composition (Section 2,
$Y=0.25$, $\ce{[Fe/H]}=-0.75$).

The upper panels of Figure 4 show the turn-off point and the subgiant branch
of the ridgeline isochrone, evaluated at $4$ distinct ages and overplotted on
the color-magnitude distribution of the archival photometry (Section 3) in the
optical (left) and near infrared (right). The isochrones were calculated using
the nominal metallicity and helium mass fraction, as well as
$E(B-V)=\qty{0.04}{mag}$, in agreement with Percival et al. (2002). In the
figure, the $\qty{11.5}{Gyr}$ isochrone is the only one that matches both the
color of the turn-off point and the luminosity of the subgiant branch within
one standard deviation. Furthermore, the presented fit firmly fixes the
reddening at $\qty{0.04}{mag}$, as lower values would result in the isochrone
being “too blue” at the turn-off point, while higher values would make the
isochrone “too faint” at the tip of the subgiant branch. This age estimate
broadly agrees with the literature (e.g. VandenBerg 2000; also see Rennó et
al. 2020 for a compilation of recent age measurements).

To explore the effect of the helium mass fraction and the overall metallicity,
we calculated two additional isochrones with $Y=0.28$ and $\ce{[Fe/H]}=-0.7$,
with all other parameters matching those of the ridgeline isochrone. The
results are plotted in the middle panels of Figure 4. While the enhanced
metallicity offers a marginal improvement of the fit at the turn-off point,
the corresponding decrease in luminosity of the subgiant branch reduces the
overall goodness of fit. On the other hand, the higher helium mass fraction
improves the fit at the subgiant branch, at the expense of mismatching the
color of the turn-off point by nearly two standard deviations in near
infrared. Based on these results, we chose to use the nominal metallicity and
helium mass fraction for all isochrone calculations henceforth.

The lower panels of the figure demonstrate how our choice of the mixing length
compares to two alternative values: $\alpha_{\mathrm{MLT}}=1.6$ and
$\alpha_{\mathrm{MLT}}=2.0$. The isochrones are shown at $\qty{11.5}{Gyr}$.
Around the turn-off point, the effect is practically indistinguishable from
that of the cluster age, emphasizing the limitations of using the color of the
turn-off point as an age diagnostic for globular clusters. Our earlier claim
that the $\qty{11.5}{Gyr}$ isochrone is the only one that fits the observed
scatter within one standard deviation is invalid if the mixing length is
allowed to vary as a free parameter. We found that a comparable goodness-of-
fit is obtained at $\qty{12.5}{Gyr}$ for $\alpha_{\mathrm{MLT}}=2.0$, or
$\qty{10.5}{Gyr}$ for $\alpha_{\mathrm{MLT}}=1.6$. Figure 4 appears to
indicate that the effect of the mixing length is most significant at the tip
of the subgiant branch. If the trend continues into the red giant branch, it
is possible that a better estimate of the cluster age may be obtained from
higher-mass stars. Since the degeneracy cannot be broken within the mass range
considered in this study, we chose to adopt the solar-calibrated
$\alpha_{\mathrm{MLT}}=1.82$ for all isochrones and the corresponding best-fit
age of $\qty{11.5}{Gyr}$.

###  5.2 Lower main sequence

Below the main sequence inflection point ($\sim\qty{0.5}{}$, Calamida et al.
2015), the abundances of individual chemical elements may noticeably impact
the shape of the CMD. In this study, we consider the variations in $11$
element abundances: [C/Fe], [N/Fe], [O/Fe], [Na/Fe], [Mg/Fe], [Al/Fe],
[Si/Fe], [S/Fe], [Ca/Fe], [Ti/Fe], and [V/Fe].

The effect of atomic abundances is twofold. First, variations in abundances
displace the chemical equilibrium of the atmosphere and change the
corresponding opacity distribution, resulting in altered photometric colors.
Second, the new atmospheric structure affects the boundary conditions of the
atmosphere-interior coupling (Section 4.2), thereby offsetting the end-point
of the evolutionary track to different $T_{\mathrm{eff}}{}$ and
$\log_{10}(g){}$. The importance of the latter effect approximately correlates
with $\delta\kappa=\left|d\kappa/d\ce{[X/Fe]}\right|$ at
$T(\tau)=T_{\mathrm{eff}}$, where [X/Fe] is the abundance of the element of
interest and $\kappa$ is the Rosseland mean opacity. For the range
$\qty{0.1}{}<M<\qty{0.5}{}$, this diagnostic is by far the largest for [O/Fe]
and [Ti/Fe] ($\sim\qty{e-7}{cm^{-1}dex^{-1}}$ at $\qty{0.2}{}$) due to the
prominent TiO and H_2O absorption bands (see Figure 2). For comparison, the
next two most important elements, [C/Fe] and [Al/Fe], have
$\delta\kappa\sim\qty{2e-8}{cm^{-1}dex^{-1}}$ and
$\sim\qty{e-8}{cm^{-1}dex^{-1}}$, respectively, at the same initial stellar
mass.

The relative importance of the two effects depends on the chosen photometric
bands. For instance, the optical main sequence spectra for fixed
$T_{\mathrm{eff}}{}$ and $\log_{10}(g){}$ are only weakly sensitive to the
oxygen abundance, since H_2O bands are mostly confined to infrared
wavelengths, while the rate of TiO production in the atmosphere is primarily
determined by [Ti/Fe]. As such, the atmosphere-interior coupling alone is
responsible for over $50\%$ of the correlation between the optical colors
(F606W – F814W) and [O/Fe] for main sequence stars. This example emphasizes
that accounting for the enhancements of individual elements in the atmosphere-
interior coupling scheme is essential for the accurate photometric
determination of the chemical composition. Since 47 Tucanae is known to have a
significant spread in [O/Fe] based on spectroscopic measurements, the
corresponding distribution of photometric colors cannot be captured with
atmosphere-interior boundary conditions based on solar or even solar-scaled
abundances. Furthermore, since both [Ti/Fe] and [O/Fe] have comparable
$\delta\kappa$, and since both are considered to be $\alpha$-elements, even
solar-scaled boundary conditions with adjustable $\alpha$-enhancement would
not be adequate.

###  5.3 Abundance variation grids

We determine the final red tail, blue tail and ridgeline isochrones for 47
Tucanae by iteratively correcting the nominal element abundances, derived in
Section 2. Each iteration begins by calculating the theoretical isochrone with
fully self-consistent chemistry as detailed in Section 4. However, since the
available datasets only extend to the end of the main sequence (Section 3), we
terminate all intermediate isochrones at $T_{\mathrm{eff}}{}=\qty{3000}{}$ to
avoid unnecessary calculations of model atmospheres.

Figure 5: Photometric sensitivity plot for the initial iteration with nominal
chemistry at $M>\qty{0.1}{}$ and for the final ridgeline isochrone at
$M<\qty{0.1}{}$. The plot shows the effect of changes in the abundances of
individual elements on the photometric colors as a function of stellar mass
(or, equivalently, $T_{\mathrm{eff}}{}$) assuming that the atmosphere-interior
coupling remains constant, i.e. for fixed $T_{\mathrm{eff}}{}$ and
$\log_{10}(g){}$. Note that the horizontal axis is categorical, i.e. the
labels are evenly spaced regardless of their values. The data shown in this
figure are available as a machine-readable table in the digital version of
this publication.

To determine the corrections in the element abundances for the next iteration,
we first compute an abundance variation grid of model atmospheres for the
current isochrone. The grid spans $5$ initial stellar masses between
$\qty{0.13}{}$ and $\qty{0.7}{}$ that sample the main sequence of the cluster
with approximately even intervals in the CMD. At this stage, we assume that
the effect of chemistry on the atmosphere-interior coupling is negligible and
compute new model atmospheres (Section 4.3) for the $T_{\mathrm{eff}}{}$ and
$\log_{10}(g){}$ of the chosen masses, based on the parameter relationships of
the current isochrone. For each stellar mass, we calculate $22$ new model
atmospheres with each of the $11$ elements considered in this study perturbed
by $\qty{+0.5}{dex}$ and $\qty{-0.5}{dex}$, one element at a time. A
photometric sensitivity plot (Figure 5) is then produced that shows the effect
of each element on the relevant photometric colors. Since the abundance
variation grids are based on the assumption of constant atmosphere-interior
coupling, we also calculate $\delta\kappa$ for each element to estimate our
confidence in the derived correlation between colors and abundances.

Figure 6: The final best fit, blue tail and red tail isochrones, overplotted
on the observed photometric spread in the optical (top) and the near infrared
(bottom) bands. The parameters of the isochrones are given in Table 1. The
isochrones shown in this figure are available as machine-readable tables in
the digital version of this publication.

The abundance corrections are then determined manually based on the current
discrepancies between the isochrones and the data, the photometric sensitivity
plot and estimated adjustments to the photometric sensitivity based on the
$\delta\kappa$ values and our experience with earlier iterations. At each
iteration, we follow the reductionist approach of adopting the smallest
possible number of corrections in the chemical composition to reproduce the
data. In principle, the derivation of abundance corrections may be automated
by introducing a quantitative goodness-of-fit criterion, similar to the one
derived in our previous work for fitting a single isochrone to the CMD
(Gerasimov et al., 2022a). However, since there is no straightforward way to
extend that criterion to simultaneous fitting of multiple isochrones, and
since the automated routine must be “trained” to take advantage of the
abundance variation grids based on experience, we chose to use human input at
every iteration instead.

A satisfactory fit was obtained after $10$ iterations. The final isochrones
were extended to $T_{\mathrm{eff}}{}\gtrsim\qty{700}{}$ to reach the
substellar regime, as required for our further analysis. For the final
ridgeline isochrone, an additional abundance variation grid was computed for
$7$ initial masses below the end of the main sequence at $\qty{0.05}{}\leq
M\leq\qty{0.08}{}$. The derived photometric sensitivity is shown in Figure 5.
We determined that corrections in [Ti/Fe] and [O/Fe] are sufficient to
reproduce the photometric spread in both optical and infrared colors, while
variations in other abundances compared to their spectroscopic means do not
offer a noticeable improvement to the fit. While the value of [Ti/Fe] needed
to be offset from its nominal value to fit the data, we found that the
observations are most consistent with a constant value of [Ti/Fe] across all
three final isochrones.

The final isochrones are plotted in Figure 6, while their best-fit properties
are summarized in Table 1. In near infrared bands, the photometric spread
noticeably overflows the red tail isochrone around the main sequence knee
(F160W $\sim\qty{19}{mag}$). A similar red excess in the optical CMD is also
seen at the corresponding evolutionary phase (F814W $\sim\qty{20.5}{mag}$),
despite the inverted [O/Fe]-color relationship. This prominent discrepancy
between the data and the model isochrones cannot be accounted for by varying
any of the considered elements.

Table 1: Best-fit parameters of the final isochrones

Parameter | Blue tail | Ridgeline | Red tail
---|---|---|---
$E(B-V)$ [$\qty{}{mag}$] | $0.04$ | $0.04$ | $0.04$
$\alpha_{\mathrm{MLT}}$ | $1.82$ | $1.82$ | $1.82$
Age [$\qty{}{Gyr}$] | $11.5$ | $11.5$ | $11.5$
[Fe/H] [$\qty{}{dex}$] | $-0.75$ | $-0.75$ | $-0.75$
$Y$ | $0.25$ | $0.25$ | $0.25$
[Ti/Fe] [$\qty{}{dex}$] | $0.64$ | $0.64$ | $0.64$
[O/Fe] [$\qty{}{dex}$] | $0.48$ | $0.33$ | $0.03$

Note. — The abundances of individual elements unlisted in the table are set to
their nominal values, derived in Section 2 and summarized in Appendix A.

##  6 Analysis

###  6.1 Oxygen abundance

Figure 6 demonstrates that the near infrared photometry of the lower main
sequence is particularly sensitive to the oxygen abundance of 47 Tucanae.
Here, we calculate the lower main sequence distribution of [O/Fe] in the
cluster. Each of the three isochrones shown in the figure defines the
relationship between the initial mass of the star and its position in the CMD
space for a given oxygen abundance ($\ce{[O/Fe]}=\qty{0.48}{dex}$ for the blue
tail isochrone, $\ce{[O/Fe]}=\qty{0.33}{dex}$ for ridgeline and
$\ce{[O/Fe]}=\qty{0.03}{dex}$ for red tail; see Table 1). Equivalent
relationships for other values of [O/Fe] may be approximated by interpolating
or extrapolating the calculated isochrones. It is therefore possible to define
a curvilinear transformation from the CMD space in Figure 6 to the
$M-\ce{[O/Fe]}$ space.

Figure 7: Left: distribution histograms of [O/Fe] in 47 Tucanae in different
initial mass bins. The flattening of the highest-mass distribution is likely
an artifact of the poor isochrone fit in the vicinity of the main sequence
knee (see Section 5.3). The lowest-mass distribution is contaminated by SMC
members. The horizontal dashed lines indicate the range of masses where the
analysis of the oxygen distribution was carried out. Right: the near infrared
CMD of 47 Tucanae, transformed into the initial mass – oxygen abundance space
using the calculated isochrones. The transformed red tail, ridgeline and blue
tail isochrones are indicated with the red, black and blue vertical lines,
respectively. Bottom: combined [O/Fe] distribution at
$\qty{0.15}{}<M<\qty{0.4}{}$, overplotted on the distribution of spectroscopic
measurements in the literature (Section 2). The measurements of mass and
oxygen abundance shown in this figure are available as a machine-readable
table in the digital version of this publication.

To obtain the desired transformation, we first generated a $1000\times 300$
regular grid of synthetic stars at $1000$ evenly spaced initial masses from
the minimum ($\approx\qty{0.045}{}$) to the maximum ($\approx\qty{0.875}{}$)
value within the range of all three isochrones, and $300$ evenly spaced values
of [O/Fe] from $\qty{-0.57}{dex}$ to $\qty{0.78}{dex}$. The chosen range of
oxygen abundances corresponds to tripled differences between the ridgeline
oxygen abundance and the other two isochrones. Each synthetic star was
projected onto the near infrared color-magnitude plane by, first, linearly
interpolating the CMD of each of the three isochrones to the mass of the
synthetic star and, second, linearly interpolating and extrapolating the
resulting three points to the oxygen abundance of the synthetic star. The
curvilinear transformation of the observed CMD of the cluster was carried out
by identifying the closest synthetic star to each real observation in the
color-magnitude space, and assigning the mass and [O/Fe] values of the closest
synthetic star to the real star. This method is unreliable above the main
sequence knee, where all three isochrones merge and the curvilinear
transformation becomes non-unique. We carry out this analysis below the knee
($\qty{0.1}{}\leq M\leq\qty{0.5}{}$), where the isochrones are well-separated,
and the transformation is accurate to $\qty{e-3}{}$ in mass and
$\qty{e-2}{dex}$ in oxygen abundance, as estimated by applying the
transformation to ${10}^{4}$ additional synthetic stars with randomly chosen
masses and oxygen abundances. We further restrict our analysis to stars that
satisfy $\qty{-0.27}{dex}<\ce{[O/Fe]}<\qty{0.63}{dex}$ (doubled differences
between the calculated isochrones) to avoid excessive extrapolation and to
discard the occasional outlying measurements that fell beyond the region of
the CMD covered with synthetic stars.

Table 2: Comparison of synthetic and real distributions of spectroscopically
inferred [O/Fe]

Statistic | Synthetic | Real
---|---|---
Mean | $0.232$ | $\pm 0.024$ | $0.174$
Standard deviation | $0.253$ | $\pm 0.015$ | $0.214$
5th percentile | $-0.198$ | $\pm 0.047$ | $-0.170$
25th percentile | $0.058$ | $\pm 0.036$ | $0.030$
75th percentile | $0.414$ | $\pm 0.029$ | $0.320$
95th percentile | $0.616$ | $\pm 0.038$ | $0.476$

Note. — All values are quoted in dex. The indicated uncertainty of the
synthetic values was calculated as the standard deviation across $10^{5}$
Monte-Carlo trials.

The result of the transformation is shown in the right panel of Figure 7. The
distribution of [O/Fe] in selected mass bins is shown in the left panel of the
figure. We note that the highest mass bin ($\qty{0.4}{}<M<\qty{0.5}{}$) is an
unreliable indicator of the oxygen distribution due to the proximity of the
main sequence knee, where the isochrone fit is noticeably poorer. The lower
half of the lowest mass bin ($\qty{0.1}{}<M<\qty{0.15}{}$) is visibly affected
by the contamination of the CMD by SMC members. The oxygen abundances within
the remaining mass range ($\qty{0.15}{}<M<\qty{0.4}{}$) were combined into a
single oxygen abundance distribution in the lower panel of Figure 7. A
histogram of the spectroscopic measurements of higher-mass stars
($M\gtrsim\qty{0.9}{}$) is plotted in the same panel for reference.

To compare the spectroscopic and photometric abundances quantitatively, we
began with a null-hypothesis that (1) the photometric distribution is free of
systematic errors and represents the true chemistry in the lower main
sequence, and (2) photometric and spectroscopic abundances are identical, i.e.
there is no dependence of [O/Fe] on stellar mass. The spectroscopic
distribution of oxygen, as described in Section 2, consists of $117$
individual measurements with published uncertainties. We generated $10^{5}$
synthetic equivalents of the spectroscopic dataset by drawing $117$ random
measurements from the photometric distribution for each trial and applying
simulated Gaussian noise according to the published spectroscopic
uncertainties. For each synthetic set, we estimated the mean, standard
deviation, 5th, 25th, 75th and 95th percentiles. The averages and spreads in
those statistics across all $10^{5}$ Monte-Carlo trials are summarized in
Table 2, alongside the equivalent values calculated for the real spectroscopic
dataset. We note that the 25th and 75th percentiles fall within the range of
[O/Fe] covered by the model isochrones and are, therefore, immune to
extrapolation errors. All other statistics listed in the table are sensitive
to the systematic offsets introduced by extrapolation.

The null-hypothesis was rejected, since some of the statistical parameters of
the real spectroscopic distribution fall outside the corresponding synthetic
ranges by over $\sim 3$ standard deviations, most notably including the
extrapolation-immune 75th percentile. If the systematic errors in the
photometrically inferred abundances are the dominant contributor to the
estimated discrepancy, their magnitude appears to vary from
$\approx\qty{0.03}{dex}$ at the oxygen-poor tail of the distribution (inferred
from the 5th percentile) to $\approx\qty{0.15}{dex}$ at the oxygen-rich tail
(inferred from the 95th percentile). We note that the latter value is
comparable to the average spectroscopic uncertainty of $\qty{0.147}{dex}$
across the $117$ measurements.

Alternatively, a genuine variation in chemistry with stellar mass may be
responsible for the discrepancy between spectroscopic and photometric values.
Under the standard assumption that the enriched population of 47 Tucanae is
oxygen-deficient (Dickens et al., 1991), the values of the higher percentiles
are expected to be primarily determined by the oxygen content of the
primordial population, while the values of the lower percentiles would be
mostly set by [O/Fe] of the enriched population. As discussed in Section 1,
the chemistry of the enriched population may depend on stellar mass in
concurrent formation models. We, however, disfavor this explanation, since the
spectroscopic/photometric discrepancy in Table 2 appears largest at the
oxygen-rich (primordial) tail of the [O/Fe] distribution, rather than the
anticipated oxygen-poor (enriched) tail.

Both photometric and spectroscopic distributions in the lower panel of Figure
7 exhibit a clear negative skewness (i.e. the oxygen-poor tail is longer than
the oxygen-rich one). To assess the statistical significance of this feature,
we calculated the Fisher-Pearson coefficient of skewness for both
distributions, obtaining $-0.21\pm 0.03$ and $-0.43\pm 0.19$ for the
photometric and spectroscopic cases respectively. Hence, spectroscopic and
photometric skewness estimates are consistent with each other, and indicate
that the distribution of [O/Fe] in 47 Tucanae is indeed negatively skewed.
While both skewness estimates have high confidence, the photometric value is
vastly more statistically significant than its spectroscopic counterpart
($8.4$ vs $\sim 2.2$ standard deviations, respectively), due to the number of
photometric measurements ($4862$) exceeding the number of spectroscopic
measurements ($117$) by over an order of magnitude.

###  6.2 Mass function

The mass functions of globular clusters retain a footprint of their dynamical
evolution (Sollima & Baumgardt, 2017; Baumgardt & Sollima, 2017) and may vary
between populations if they have distinct kinematic properties. Furthermore,
an estimate of the mass function is required to predict the luminosity
function and the CMD of the substellar sequence. In this sub-section, we seek
a mass function that can reproduce the observed F160W magnitude distribution
of 47 Tucanae, assuming that the isochrones calculated earlier and the oxygen
abundance distribution derived above accurately capture the properties of the
cluster. The choice of the photometric band is motivated by the fact that the
oxygen abundances were inferred from the near infrared CMD and the fact that
lower main sequence members are brighter in F160W than in F110W (Figure 6). We
also assumed that the mass function ($\xi$) is of the form of a broken power
law (Kroupa, 2001):

$\xi(M)\propto\begin{cases}M^{-\beta},&\text{if }M>M_{\mathrm{bp}}\\\
M^{-\gamma},&\text{if }M\leq M_{\mathrm{bp}}\\\ \end{cases}$ (2)

Here, $\beta$ and $\gamma$ are the power law indices in the high- and low-mass
regimes respectively, separated by the break-point stellar mass
$M_{\mathrm{bp}}$. We keep all three parameters as free variables. The
magnitude function, $\phi(m)$, where $m$ is the apparent magnitude of the star
(in our case, F160W), is related to the mass function as in Equation 3:

$\phi(m)=-\xi\left(M(m)\right)\frac{dM(m)}{dm}$ (3)

In Equation 3, we treat the stellar mass, $M$, as a function of magnitude,
$m$, as given by the mass-magnitude relationship provided by the isochrone. In
practice, the mass-magnitude relationship is sampled at a set of masses,
$M_{i}$, and magnitudes, $m_{i}$, corresponding to the initial masses of the
calculated evolutionary models. We therefore rewrite Equation 3 in the form of
finite differences:

$\phi\left(\frac{m_{i}+m_{i+1}}{2}\right)=-\xi\left(\frac{M_{i}+M_{i+1}}{2}\right)\frac{M_{i+1}-M_{i}}{m_{i+1}-m_{i}}$
(4)

Figure 8: The observed magnitude function of 47 Tucanae in F160W with the
best-fit theoretical magnitude function in magenta. Additional magnitude
functions that use only one of the theoretical isochrones instead of the
abundance-weighted average are shown for comparison. The displayed error bars
are Poisson counting errors. The observed counts in the four faintest bins are
excluded from this analysis due to contamination by SMC members.

Here, the magnitude function is evaluated at the mid-points between the
adjacent evolutionary models on the mass grid. We calculated three distinct
magnitude functions for each of the three final isochrones (red tail, blue
tail and ridgeline). The magnitude function for an arbitrary [O/Fe] can be
obtained by linearly interpolating/extrapolating between the three calculated
magnitude functions to the desired oxygen abundance. The combined magnitude
function for the entire cluster was computed as the average of individual
magnitude functions, evaluated for the oxygen abundance of every star that
satisfies $\qty{0.15}{}<M<\qty{0.4}{}$ and
$\qty{-0.27}{dex}<\ce{[O/Fe]}<\qty{0.63}{dex}$ (the same restrictions as in
the lower panel of Figure 7). Finally, the magnitude function was binned into
$20$ uniform bins in the range $17<m<24$ using trapezoid integration. The four
faintest bins with bin centers at F160W $>22.5$ were excluded from our
analysis due to contamination by SMC members and potential incompleteness of
the photometric sample.

To determine the parameters of the mass function, $\beta$, $\gamma$ and
$M_{\mathrm{bp}}$, we applied the same binning to the observed magnitude
function in F160W and fitted the theoretical magnitude function to the result,
using least squares regression, weighted by the Poisson errors in each bin.
The observed magnitude function and the theoretical best fit are shown in
Figure 8. The best-fit parameters of the mass function are $\beta=1.75\pm
0.15$, $\gamma=0.08\pm 0.04$ and $M_{\mathrm{bp}}=\qty{0.47\pm 0.02}{}$. In
addition to the combined magnitude function, three more theoretical magnitude
functions were calculated for the same parameters as the best fit, but using
only one of the three isochrones instead of a weighted average of all three.
As seen in Figure 8, the scatter among the theoretical magnitude functions is
consistent with the Poisson counting errors in the corresponding magnitude
bins, suggesting that for main sequence members the distribution of oxygen
abundance cannot be inferred from the observed luminosity function of the
cluster.

###  6.3 Brown dwarfs in 47 Tucanae

Figure 9: Predicted magnitude distribution of stars and brown dwarfs in 47
Tucanae at three different ages. The bright and faint peaks represent the main
sequence and the brown dwarfs of the cluster, respectively, with a clear
stellar/substellar gap in between. The distribution is normalized to the
number of stars in the dataset (Section 3). An animated version of this figure
is available in the digital version of the publication. The animation lasts
$6$ seconds and includes $30$ frames. The animation shows the evolution of the
magnitude function from $\qty{0.1}{Gyr}$ to $\qty{13.5}{Gyr}$. Initially, the
magnitude function only has one peak around $m\sim 20$ corresponding to the
main sequence of the cluster. The stellar/substellar gap forms at $m\sim 23.5$
around $\qty{0.5}{Gyr}$. Over time, it gradually deepens, widens, and shifts
to fainter magnitudes. Figure 10: Predicted CMD of 47 Tucanae from the
subgiant branch to the substellar sequence. Isomass lines are shown for
selected masses. The color scatter in the figure is derived from the inferred
distribution of [O/Fe] near the end of the main sequence. The inset sub-figure
in the upper right corner shows the mass-magnitude relationship for the
ridgeline isochrone with the hydrogen-burning limit highlighted. The
isochrones used to construct the CMD in this figure are available as machine-
readable tables in the digital version of this publication.

In this sub-section, we predict the colors, magnitudes and number densities of
the brown dwarfs in 47 Tucanae, as they may be observed with JWST NIRCam. The
JWST color-magnitude space considered in this study is F322W2 vs F150W2 –
F322W2, since brown dwarfs are most likely to be detected in wide bands at
long wavelengths due to their faint magnitudes and red colors. Our predictions
are based on the assumption that the mass function (Equation 2 with the best-
fit parameters in Figure 8), the calculated isochrones (Figure 6) and the
inferred distribution of chemical abundances (Figure 7) remain unchanged in
the substellar regime.

We re-calculated the synthetic magnitude function of 47 Tucanae in $30$ evenly
spaced magnitude bins within range of all three theoretical isochrones, using
the best-fit mass function parameters. The calculations were carried out for a
variety of ages between $\qty{0.1}{Gyr}$ and $\qty{13.5}{Gyr}$, including
$\qty{11.5}{Gyr}$ (best-fit age). The requirement to vary the age of the
cluster for this part of the analysis necessitated the use of interpolated
bolometric corrections for synthetic photometry instead of atmosphere
hammering. The resulting magnitude functions for three of the considered ages
are shown in the static preview of Figure 9. At those ages, the magnitude
function exhibits a clear stellar/substellar gap at $m\sim 26$, where the
brown dwarf density per magnitude is reduced by nearly an order of magnitude
(compared to the maximum density within the modelled range that is predicted
to occur around $m=27.5$). While the lower main sequence is virtually
unaffected by the cluster age, both the turn-off point and the substellar
sequence undergo noticeable evolution. The equivalent magnitude functions for
all other ages are available as an animated figure in the digital version of
this paper.

A synthetic CMD based on the inferred oxygen abundance in Figure 7, the
isochrones in Figure 6 and the predicted mass function at $\qty{11.5}{Gyr}$ is
shown in Figure 10. Lines of constant mass (isomass lines) for selected masses
and $\qty{-0.27}{dex}<\ce{[O/Fe]}<\qty{0.63}{dex}$ are also shown for
reference. Since the isochrones of the substellar sequence in 47 Tucanae
appear almost exactly parallel to the isomass lines, the effect of evolution
on the brown dwarf colors is highly degenerate with [O/Fe]. This presents a
potential challenge to the use of brown dwarfs as chemical tracers, since the
masses of individual brown dwarfs (and, hence, their evolutionary phases) are
not a priori known. The degeneracy also reduces the observed photometric
scatter, making the substellar sequence far narrower than the main sequence in
the CMD space.

The predicted mass-magnitude relationship for the ridgeline isochrone is shown
in an inset sub-figure of Figure 10 with the hydrogen-burning limit (HBL)
highlighted. We define the HBL as the mass, at which the energy output from
nuclear reactions contributes $50\%$ of the stellar luminosity at
$\qty{11.5}{Gyr}$. In this case, the HBL was computed as $\qty{0.074}{}$.

##  7 Discussion and conclusion

We developed a general method for determining the chemical compositions, ages
and mass functions of globular clusters from the observed CMDs in optical and
near infrared bands. Our method relies on state-of-the-art model atmospheres
and evolutionary models that are fully self-consistent and incorporate the
full set of non-solar abundances in every component, including the interior
structure and evolution, nuclear processing, atmosphere-interior coupling,
atmospheric structures and spectral synthesis. Our modelling framework was
applied to the brightest mono-metallic globular cluster 47 Tucanae. We
reproduce for the first time the observed scatter in the photometric colors of
main sequence stars without any a priori assumptions of its magnitude. We also
provide the first measurements of chemical compositions for individual stars
of the lower main sequence in a globular cluster (Figure 7). An extension of
our models to the substellar regime predicts the expected colors and
magnitudes of brown dwarfs in the cluster, reproducing the anticipated
stellar/substellar gap. The predicted brown dwarf CMD for the first time
incorporates the inferred distribution of chemical abundances. The best-fit
parameters of 47 Tucanae determined in this study are listed in Table 1. The
key findings are as follows:

  * •

The photometric scatter of 47 Tucanae in both optical and near infrared bands
can be reproduced with unaltered average spectroscopic abundances for all
elements (Appendix A) with the exception of [O/Fe] and [Ti/Fe]. Constraining
the abundances of elements other than oxygen and titanium from the CMDs
considered in this study is challenging due to the subdominant effect of those
elements on stellar atmospheres, as illustrated in Figure 5.

  * •

The observed photometric scatter is predominantly driven by [O/Fe]. Our best-
fit model CMD is consistent with the absence of star-to-star variations in
[Ti/Fe], in agreement with the spectroscopic measurements and the mono-
metallic nature of 47 Tucanae.

  * •

The photometric distribution of the oxygen abundance estimated in this study
(Figure 7) at $M<\qty{0.4}{}$ is statistically consistent with the
spectroscopic distribution, inferred from evolved stars, under the assumption
of systematic errors in the photometric estimates of order $\qty{0.15}{dex}$.
These errors are comparable to the experimental uncertainties in available
spectroscopic measurements, demonstrating that lower main sequence CMDs of
globular clusters may be used in chemical analyses.

  * •

The discrepancy between the spectroscopic and photometric [O/Fe] values
appears largest at the oxygen-rich (primordial) tail of the distribution. If a
genuine mass dependence of chemistry makes a significant contribution to the
observed discrepancy in addition to the suspected systematic offset, the
variation of [O/Fe] with stellar mass would be most prominent in the enriched
population of the cluster instead of the primordial population. Since this
conclusion contradicts the theoretical expectation, we disfavor this
explanation and suggest that the mass dependence of chemistry makes a
subdominant contribution compared to the systematic offset, if any. We
therefore expect the difference in [O/Fe] between the lower main sequence and
evolved stars to be less than the estimated systematic offset of
$~{}\qty{0.15}{dex}$.

  * •

Corroborating our previous study of $\omega$ Centauri, we confirm that the
JWST CMD of 47 Tucanae is expected to have a stellar/substellar gap below the
end of the main sequence, followed by a large number of brown dwarfs (Figure
9). The gap occurs around F322W2$\sim\qty{26}{mag}$ and the maximum density of
brown dwarfs within the modelling range is attained at
F322W2$\sim\qty{27.5}{mag}$. Our evolutionary models suggest that
stellar/substellar gaps form around the age of $\qty{0.5}{Gyr}$ in globular
clusters similar to 47 Tucanae, and get deeper and wider over time.

  * •

While our models indicate that the colors of brown dwarfs are highly sensitive
to chemical abundances (Figure 5), the apparent degeneracy between their
isomass lines and isochrones (Figure 10) will make the inference of oxygen
abundance challenging in the considered color-magnitude space (F322W2 vs
F150W2 – F322W2). On the other hand, the brown dwarf magnitude function will
be sensitive to chemistry, introducing additional systematic biases to age
estimates inferred from the stellar/substellar gap.

The constraints on the mass dependence of the oxygen abundance inferred in
this study restrict the allowed parameter space of concurrent formation MP
models, in agreement with previous studies of the photometric spread near the
end of the main sequence (Ziliotto et al., 2023; Milone, 2023). However, a
detailed model of pollution as a function of stellar mass is required to carry
out a thorough evaluation of proposed concurrent formation models, which may
be challenging to derive due to the uncertain properties and evolution of
circumstellar disks.

Future observations of the substellar sequences in nearby globular clusters,
including 47 Tucanae, will extend the mass baseline of photometric analysis by
nearly an order of magnitude and explore the abundances of other elements such
as carbon and nitrogen that are far more important in the atmospheres of brown
dwarfs than main sequence stars (Figure 5). The aforementioned isomass-
isochrone degeneracy suggests that the luminosity function of brown dwarfs in
globular clusters needs to be considered in conjunction with the CMD when
deriving photometric chemical abundances. We further emphasize that JWST and
HST bands other than the ones considered here may be more sensitive to
chemical variations in substellar atmospheres, necessitating a follow-up study
that explores all possible color combinations rather than a small subset of
filters used in this paper.

We note that the predicted substellar CMD of 47 Tucanae in Figure 10 was
calculated under the assumption that the variation in [O/Fe] alone is
sufficient to reproduce the photometric scatter. While this assumption was
shown to be accurate for the main sequence stars, it likely has limited
validity in the substellar regime and a more accurate theoretical CMD may be
obtained by calculating additional isochrones that incorporate the
spectroscopically inferred scatter in other elements.

If the spectroscopic/photometric discrepancy discussed here and in Section 6.1
is primarily caused by systematic offsets, their nature must be investigated
in a future study. In particular, the accuracy of interpolation and
extrapolation of best-fit theoretical isochrones can be verified by a denser
sampling of [O/Fe] with calculated stellar models. The accuracy of the models
themselves can be assessed by comparing synthetic spectra to spectroscopic
observations of nearby metal-poor stars and brown dwarfs.

We acknowledge the funding support from Hubble Space Telescope (HST) Program
GO-15096, provided by NASA through a grant from the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Incorporated, under NASA contract NAS5-26555. The computational
demand of this study was met by the Extreme Science and Engineering Discovery
Environment (XSEDE), supported by NSF grant ACI-1548562. Some of the software
used in this study was produced with assistance from ChatGPT 3.5. ChatGPT is a
large language model developed and maintained by OpenAI. This work was
conducted at University of California San Diego, which was built on the
unceded territory of the Kumeyaay Nation. Today, the Kumeyaay people continue
to maintain their political sovereignty and cultural traditions as vital
members of the San Diego community. We acknowledge their tremendous
contributions to our region and thank them for their stewardship.

##  Appendix A Nominal composition

Table 3 lists the nominal composition of 47 Tucanae adopted in this work,
based on spectroscopic observations of giant and sub-giant members of the
cluster in the literature. The derivation of the abundances listed in the
table is described in Section 2. The physical member-to-member spread provided
in the table ($s^{(X)}$) was calculated using Equation 1.

Table 3: Average chemical composition of 47 Tucanae inferred from
spectroscopic measurements in literature

| Best estimate | Uncertainty | Physical spread | # of measurements | Reference(s)
---|---|---|---|---|---
[Fe/H] | $-0.75$ | $0.01$ | $-$ | $281$ | $(1)(2)(3)$
[Al/Fe] | $0.35$ | $0.01$ | $0.06$ | $195$ | $(1)(2)(3)$
[Ba/Fe] | $0.25$ | $0.07$ | $-$ | $13$ | $(1)$
[C /Fe] | $-0.25$ | $0.01$ | $0.10$ | $70$ | $(3)$
[Ca/Fe] | $0.26$ | $0.01$ | $-$ | $162$ | $(1)(2)$
[Ce/Fe] | $-0.04$ | $0.10$ | $-$ | $11$ | $(1)$
[Co/Fe] | $-0.00$ | $0.02$ | $-$ | $13$ | $(1)$
[Cr/Fe] | $-0.04$ | $0.03$ | $-$ | $13$ | $(1)$
[Cu/Fe] | $-0.14$ | $0.10$ | $-$ | $13$ | $(1)$
[Dy/Fe] | $0.70$ | $0.07$ | $-$ | $13$ | $(1)$
[Eu/Fe] | $0.44$ | $0.01$ | $-$ | $150$ | $(1)(2)$
[La/Fe] | $0.20$ | $0.01$ | $-$ | $144$ | $(1)(2)$
[Mg/Fe] | $0.35$ | $0.01$ | $-$ | $87$ | $(1)(3)$
[Mn/Fe] | $-0.19$ | $0.04$ | $-$ | $13$ | $(1)$
[Mo/Fe] | $0.55$ | $0.04$ | $-$ | $13$ | $(1)$
[N /Fe] | $0.85$ | $0.05$ | $0.33$ | $54$ | $(3)$
[Na/Fe] | $0.27$ | $0.01$ | $0.15$ | $236$ | $(1)(2)(3)$
[Nd/Fe] | $0.04$ | $0.07$ | $-$ | $13$ | $(1)$
[Ni/Fe] | $-0.06$ | $0.01$ | $-$ | $174$ | $(1)(2)(3)$
[O /Fe] | $0.18$ | $0.02$ | $0.16$ | $117$ | $(1)(2)$
[Pr/Fe] | $-0.03$ | $0.06$ | $-$ | $13$ | $(1)$
[Ru/Fe] | $0.50$ | $0.04$ | $-$ | $13$ | $(1)$
[Sc/Fe] | $0.21$ | $0.05$ | $-$ | $13$ | $(1)$
[Si/Fe] | $0.32$ | $0.01$ | $-$ | $222$ | $(1)(2)(3)$
[Ti/Fe] | $0.34$ | $0.01$ | $-$ | $163$ | $(1)(2)$
[V /Fe] | $0.17$ | $0.04$ | $-$ | $13$ | $(1)$
[Y /Fe] | $0.07$ | $0.05$ | $-$ | $13$ | $(1)$
[Zn/Fe] | $0.26$ | $0.04$ | $-$ | $13$ | $(1)$
[Zr/Fe] | $0.41$ | $0.06$ | $0.10$ | $13$ | $(1)$

Note. — All values are quoted in dex with respect to solar abundances. (1) –
Thygesen et al. (2014). (2) – Cordero et al. (2014). (3) – Marino et al.
(2016).

## References

  * Alabi et al. (2016) Alabi, A. B., Forbes, D. A., Romanowsky, A. J., et al. 2016, MNRAS, 460, 3838, doi: 10.1093/mnras/stw1213
  * Alabi et al. (2017) —. 2017, MNRAS, 468, 3949, doi: 10.1093/mnras/stx678
  * Allard et al. (2001) Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357, doi: 10.1086/321547
  * Allard et al. (2011) Allard, F., Homeier, D., & Freytag, B. 2011, in Cool Stars 16, Vol. 448, 91, doi: 10.48550/arXiv.1011.5405
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765, doi: 10.1098/rsta.2011.0269
  * Allard et al. (2014) Allard, F., Homeier, D., & Sharp, C. 2014, Star, Brown Dwarf and Planet Simulator, ENS de Lyon.  https://phoenix.ens-lyon.fr/simulator/index.faces
  * Anderson (1997) Anderson, A. J. 1997, PhD thesis, UNIVERSITY OF CALIFORNIA, BERKELEY
  * Anderson (2002) Anderson, J. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 265, Omega Centauri, A Unique Window into Astrophysics, ed. F. van Leeuwen, J. D. Hughes, & G. Piotto, 87
  * Anderson & van der Marel (2010) Anderson, J., & van der Marel, R. P. 2010, ApJ, 710, 1032, doi: 10.1088/0004-637X/710/2/1032
  * Andrievsky et al. (2008) Andrievsky, S. M., Spite, M., Korotin, S. A., et al. 2008, A&A, 481, 481, doi: 10.1051/0004-6361:20078837
  * Angelou et al. (2011) Angelou, G. C., Church, R. P., Stancliffe, R. J., Lattanzio, J. C., & Smith, G. H. 2011, ApJ, 728, 79, doi: 10.1088/0004-637X/728/2/79
  * Arkelyan & Pilipenko (2022) Arkelyan, N. R., & Pilipenko, S. V. 2022, Astronomy Reports, 66, 191, doi: 10.1134/S1063772922030015
  * Baade (1944) Baade, W. 1944, ApJ, 100, 137, doi: 10.1086/144650
  * Bailin (2019) Bailin, J. 2019, ApJS, 245, 5, doi: 10.3847/1538-4365/ab4812
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, Astronomy & Astrophysics, 577, A42, doi: 10.1051/0004-6361/201425481
  * Barber et al. (2013) Barber, C. B., Dobkin, D. P., & Huhdanpaa, H. 2013, Qhull: Quickhull algorithm for computing the convex hull, Astrophysics Source Code Library, record ascl:1304.016.  http://ascl.net/1304.016
  * Barber et al. (2006) Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087, doi: 10.1111/j.1365-2966.2006.10184.x
  * Bastian et al. (2013) Bastian, N., Lamers, H. J. G. L. M., de Mink, S. E., et al. 2013, MNRAS, 436, 2398, doi: 10.1093/mnras/stt1745
  * Bastian & Lardo (2015) Bastian, N., & Lardo, C. 2015, MNRAS, 453, 357, doi: 10.1093/mnras/stv1661
  * Bastian & Lardo (2018) —. 2018, ARA&A, 56, 83, doi: 10.1146/annurev-astro-081817-051839
  * Baum (1952) Baum, W. A. 1952, AJ, 57, 222, doi: 10.1086/106758
  * Baumgardt & Sollima (2017) Baumgardt, H., & Sollima, S. 2017, MNRAS, 472, 744, doi: 10.1093/mnras/stx2036
  * Bedin et al. (2004) Bedin, L. R., Piotto, G., Anderson, J., et al. 2004, ApJ, 605, L125, doi: 10.1086/420847
  * Bedin et al. (2021) Bedin, L. R., Anderson, J., Apai, D., et al. 2021, The faintest and coolest stars in the two closest globulars, JWST Proposal. Cycle 1
  * Bergemann (2011) Bergemann, M. 2011, MNRAS, 413, 2184, doi: 10.1111/j.1365-2966.2011.18295.x
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, Zeitschrift für Astrophysik, 46, 108
  * Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195, doi: 10.1093/mnras/112.2.195
  * Briley et al. (1996) Briley, M. M., Smith, V. V., Suntzeff, N. B., et al. 1996, Nature, 383, 604, doi: 10.1038/383604a0
  * Brown (2005) Brown, L. R. 2005, J. Quant. Spec. Radiat. Transf., 96, 251, doi: 10.1016/j.jqsrt.2004.12.037
  * Brown et al. (2021) Brown, S. T., Buitrago, P., Hanna, E., et al. 2021, in Practice and Experience in Advanced Research Computing (ACM), doi: 10.1145/3437359.3465593
  * Brugger (1969) Brugger, R. M. 1969, The American Statistician, 23, 32–32, doi: 10.1080/00031305.1969.10481865
  * Burgasser (2004) Burgasser, A. J. 2004, ApJS, 155, 191, doi: 10.1086/424386
  * Burke et al. (2004) Burke, C. J., Pinsonneault, M. H., & Sills, A. 2004, ApJ, 604, 272, doi: 10.1086/381242
  * Burrows et al. (1989) Burrows, A., Hubbard, W. B., & Lunine, J. I. 1989, ApJ, 345, 939, doi: 10.1086/167964
  * Caiazzo et al. (2017) Caiazzo, I., Heyl, J. S., Richer, H., & Kalirai, J. 2017, arXiv e-prints, arXiv:1702.00091.  https://arxiv.org/abs/1702.00091
  * Caiazzo et al. (2019) Caiazzo, I., Burgasser, A., Rees, J. M., et al. 2019, BAAS, 51, 521.  https://arxiv.org/abs/1903.06769
  * Caiazzo et al. (2021) Caiazzo, I., Richer, H. B., Boley, A. C., et al. 2021, Brown Dwarfs, White Dwarfs and Planetary Disks in an Ancient Stellar System, JWST Proposal. Cycle 1
  * Calamida et al. (2015) Calamida, A., Sahu, K. C., Casertano, S., et al. 2015, ApJ, 810, 8, doi: 10.1088/0004-637X/810/1/8
  * Carbon (1984) Carbon, D. F. 1984, Line blanketing. (Cambridge University Press), 395–426
  * Carpenter et al. (2006) Carpenter, J. M., Mamajek, E. E., Hillenbrand, L. A., & Meyer, M. R. 2006, ApJ, 651, L49, doi: 10.1086/509121
  * Carretta et al. (2010) Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2010, A&A, 516, A55, doi: 10.1051/0004-6361/200913451
  * Carretta et al. (2000) Carretta, E., Gratton, R. G., Clementini, G., & Fusi Pecci, F. 2000, ApJ, 533, 215, doi: 10.1086/308629
  * Castelli (2005) Castelli, F. 2005, Memorie della Societa Astronomica Italiana Supplementi, 8, 34
  * Castelli & Kurucz (2003) Castelli, F., & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, Vol. 210, A20.  https://arxiv.org/abs/astro-ph/0405087
  * Chabrier & Baraffe (1997) Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039, doi: 10.48550/arXiv.astro-ph/9704118
  * Chabrier & Baraffe (2000) —. 2000, ARA&A, 38, 337, doi: 10.1146/annurev.astro.38.1.337
  * Chen et al. (2018) Chen, S., Richer, H., Caiazzo, I., & Heyl, J. 2018, ApJ, 867, 132, doi: 10.3847/1538-4357/aae089
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, The Astrophysical Journal, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Cohen (1978) Cohen, J. G. 1978, ApJ, 223, 487, doi: 10.1086/156284
  * Cordero et al. (2014) Cordero, M. J., Pilachowski, C. A., Johnson, C. I., et al. 2014, ApJ, 780, 94, doi: 10.1088/0004-637X/780/1/94
  * Cottrell & Da Costa (1981) Cottrell, P. L., & Da Costa, G. S. 1981, ApJ, 245, L79, doi: 10.1086/183527
  * Crocker et al. (1988) Crocker, D. A., Rood, R. T., & O’Connell, R. W. 1988, ApJ, 332, 236, doi: 10.1086/166648
  * Da Costa & Demarque (1982) Da Costa, G. S., & Demarque, P. 1982, ApJ, 259, 193, doi: 10.1086/160159
  * de Lacaille (1755) de Lacaille, N. 1755, Mémoires de l’Académie Royale des Sciences, 194.  http://www.messier.seds.org/xtra/history/lac1755.html
  * de Silva et al. (2009) de Silva, G. M., Gibson, B. K., Lattanzio, J., & Asplund, M. 2009, A&A, 500, L25, doi: 10.1051/0004-6361/200912279
  * Decressin et al. (2007) Decressin, T., Charbonnel, C., & Meynet, G. 2007, A&A, 475, 859, doi: 10.1051/0004-6361:20078425
  * Denissenkov & Hartwick (2014) Denissenkov, P. A., & Hartwick, F. D. A. 2014, MNRAS, 437, L21, doi: 10.1093/mnrasl/slt133
  * D’Ercole et al. (2016) D’Ercole, A., D’Antona, F., & Vesperini, E. 2016, MNRAS, 461, 4088, doi: 10.1093/mnras/stw1583
  * D’Ercole et al. (2008) D’Ercole, A., Vesperini, E., D’Antona, F., McMillan, S. L. W., & Recchi, S. 2008, MNRAS, 391, 825, doi: 10.1111/j.1365-2966.2008.13915.x
  * Dickens & Bell (1976) Dickens, R. J., & Bell, R. A. 1976, ApJ, 207, 506, doi: 10.1086/154518
  * Dickens et al. (1991) Dickens, R. J., Croke, B. F. W., Cannon, R. D., & Bell, R. A. 1991, Nature, 351, 212, doi: 10.1038/351212a0
  * Dieball et al. (2016) Dieball, A., Allard, F., Bedin, L. R., et al. 2016, Hunting for Brown Dwarfs in Globular Clusters: Second Epoch Deep IR observations of the Globular Clusters M4, HST Proposal
  * Dieball et al. (2019) Dieball, A., Bedin, L. R., Knigge, C., et al. 2019, MNRAS, 486, 2254, doi: 10.1093/mnras/stz996
  * Doppel et al. (2021) Doppel, J. E., Sales, L. V., Navarro, J. F., et al. 2021, MNRAS, 502, 1661, doi: 10.1093/mnras/staa3915
  * Dorman et al. (1989) Dorman, B., Nelson, L. A., & Chau, W. Y. 1989, ApJ, 342, 1003, doi: 10.1086/167658
  * Dotter (2013) Dotter, A. 2013, Mem. Soc. Astron. Italiana, 84, 97, doi: 10.48550/arXiv.1307.5589
  * Dotter et al. (2015) Dotter, A., Ferguson, J. W., Conroy, C., et al. 2015, MNRAS, 446, 1641, doi: 10.1093/mnras/stu2170
  * Fitzpatrick & Massa (2007) Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320, doi: 10.1086/518158
  * Freeman & Bland-Hawthorn (2002) Freeman, K., & Bland-Hawthorn, J. 2002, ARA&A, 40, 487, doi: 10.1146/annurev.astro.40.060401.093840
  * Freeman & Rodgers (1975) Freeman, K. C., & Rodgers, A. W. 1975, ApJ, 201, L71, doi: 10.1086/181945
  * Fried & Zietz (1973) Fried, J., & Zietz, S. 1973, Physics in Medicine and Biology, 18, 550, doi: 10.1088/0031-9155/18/4/306
  * Gamow & Keller (1945) Gamow, G., & Keller, G. 1945, Reviews of Modern Physics, 17, 125, doi: 10.1103/RevModPhys.17.125
  * Gerasimov et al. (2022a) Gerasimov, R., Burgasser, A. J., Homeier, D., et al. 2022a, ApJ, 930, 24, doi: 10.3847/1538-4357/ac61e5
  * Gerasimov et al. (2022b) Gerasimov, R., Burgasser, A. J., Homeier, D., et al. 2022b, in Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 120, doi: 10.5281/zenodo.7535820
  * Gerasimov et al. (2020) Gerasimov, R., Homeier, D., Burgasser, A., & Bedin, L. R. 2020, Research Notes of the American Astronomical Society, 4, 214, doi: 10.3847/2515-5172/abcf2c
  * Gieles et al. (2018) Gieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS, 478, 2461, doi: 10.1093/mnras/sty1059
  * Gratton et al. (2019) Gratton, R., Bragaglia, A., Carretta, E., et al. 2019, A&A Rev., 27, 8, doi: 10.1007/s00159-019-0119-3
  * Gratton et al. (2015) Gratton, R., Lucatello, S., Carretta, E., & Bragaglia, A. 2015, Highlights of Astronomy, 16, 230, doi: 10.1017/S1743921314005547
  * Gratton et al. (2004) Gratton, R., Sneden, C., & Carretta, E. 2004, ARA&A, 42, 385, doi: 10.1146/annurev.astro.42.053102.133945
  * Gratton (1982) Gratton, R. G. 1982, A&A, 115, 336
  * Gratton et al. (2000) Gratton, R. G., Sneden, C., Carretta, E., & Bragaglia, A. 2000, A&A, 354, 169
  * Gratton et al. (2001) Gratton, R. G., Bonifacio, P., Bragaglia, A., et al. 2001, A&A, 369, 87, doi: 10.1051/0004-6361:20010144
  * Greenstein (1939) Greenstein, J. L. 1939, ApJ, 90, 387, doi: 10.1086/144115
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951, doi: 10.1051/0004-6361:200809724
  * Hachenberg (1939) Hachenberg, O. 1939, ZAp, 18, 49
  * Han et al. (2009) Han, S.-I., Lee, Y.-W., Joo, S.-J., et al. 2009, ApJ, 707, L190, doi: 10.1088/0004-637X/707/2/L190
  * Harding (1962) Harding, G. A. 1962, The Observatory, 82, 205
  * Harris (1974) Harris, W. E. 1974, ApJ, 192, L161, doi: 10.1086/181616
  * Harris (1991) —. 1991, ARA&A, 29, 543, doi: 10.1146/annurev.aa.29.090191.002551
  * Harris (1996) —. 1996, AJ, 112, 1487, doi: 10.1086/118116
  * Harris & Racine (1979) Harris, W. E., & Racine, R. 1979, ARA&A, 17, 241, doi: 10.1146/annurev.aa.17.090179.001325
  * Hauschildt et al. (1999) Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377, doi: 10.1086/306745
  * Hauschildt et al. (1997) Hauschildt, P. H., Baron, E., & Allard, F. 1997, ApJ, 483, 390, doi: 10.1086/304233
  * Helling et al. (2008) Helling, C., Ackerman, A., Allard, F., et al. 2008, MNRAS, 391, 1854, doi: 10.1111/j.1365-2966.2008.13991.x
  * Hesser & Bell (1980) Hesser, J. E., & Bell, R. A. 1980, ApJ, 238, L149, doi: 10.1086/183276
  * Hilker et al. (2020) Hilker, M., Baumgardt, H., Sollima, A., & Bellini, A. 2020, in Star Clusters: From the Milky Way to the Early Universe, ed. A. Bragaglia, M. Davies, A. Sills, & E. Vesperini, Vol. 351, 451–454, doi: 10.1017/S1743921319006823
  * Homeier et al. (2003) Homeier, D., Hauschildt, P. H., & Allard, F. 2003, in Brown Dwarfs, ed. E. Martín, Vol. 211, 419.  https://arxiv.org/abs/astro-ph/0209171
  * Hubble (1932) Hubble, E. 1932, ApJ, 76, 44, doi: 10.1086/143397
  * Johnson & Pilachowski (2010) Johnson, C. I., & Pilachowski, C. A. 2010, ApJ, 722, 1373, doi: 10.1088/0004-637X/722/2/1373
  * Johnson (1994) Johnson, H. R. 1994, in Molecules in the Stellar Environment, ed. U. G. Jørgensen, Vol. 428 (Springer Berlin Heidelberg), 234–249, doi: 10.1007/3-540-57747-5_47
  * Joyce & Chaboyer (2018) Joyce, M., & Chaboyer, B. 2018, ApJ, 856, 10, doi: 10.3847/1538-4357/aab200
  * Kalirai et al. (2012) Kalirai, J. S., Richer, H. B., Anderson, J., et al. 2012, AJ, 143, 11, doi: 10.1088/0004-6256/143/1/11
  * Keenan & Keller (1953) Keenan, P. C., & Keller, G. 1953, ApJ, 117, 241, doi: 10.1086/145687
  * Kraft (1994) Kraft, R. P. 1994, PASP, 106, 553, doi: 10.1086/133416
  * Krauss & Chaboyer (2003) Krauss, L. M., & Chaboyer, B. 2003, Science, 299, 65, doi: 10.1126/science.1075631
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Kurucz (1970) Kurucz, R. L. 1970, SAO Special Report, 309
  * Kurucz (2005) —. 2005, Memorie della Societa Astronomica Italiana Supplementi, 8, 14
  * Kurucz (2014) —. 2014, Model Atmosphere Codes: ATLAS12 and ATLAS9 (Springer International Publishing), 39–51, doi: 10.1007/978-3-319-06956-2_4
  * Kurucz & Avrett (1981) Kurucz, R. L., & Avrett, E. H. 1981, SAO Special Report, 391
  * Kurucz et al. (1974) Kurucz, R. L., Peytremann, E., & Avrett, E. H. 1974, Blanketed model atmospheres for early-type stars (Smithsonian Astrophysical Observatory)
  * Larkin et al. (2023) Larkin, M. M., Gerasimov, R., & Burgasser, A. J. 2023, AJ, 165, 2, doi: 10.3847/1538-3881/ac9b43
  * Lee (2012) Lee, M. G. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 458, Galactic Archaeology: Near-Field Cosmology and the Formation of the Milky Way, ed. W. Aoki, M. Ishigaki, T. Suda, T. Tsujimoto, & N. Arimoto, 291
  * Lee (1977) Lee, S. W. 1977, A&AS, 27, 381
  * Lide (2004) Lide, D. 2004, CRC Handbook of Chemistry and Physics, 85th Edition, CRC Handbook of Chemistry and Physics, 85th Ed No. v. 85 (Taylor & Francis).  https://books.google.com/books?id=WDll8hA006AC
  * Lind et al. (2011) Lind, K., Asplund, M., Barklem, P. S., & Belyaev, A. K. 2011, A&A, 528, A103, doi: 10.1051/0004-6361/201016095
  * Lind et al. (2015) Lind, K., Koposov, S. E., Battistini, C., et al. 2015, A&A, 575, L12, doi: 10.1051/0004-6361/201425554
  * Lunine et al. (1986) Lunine, J. I., Hubbard, W. B., & Marley, M. S. 1986, ApJ, 310, 238, doi: 10.1086/164678
  * MacLean et al. (2015) MacLean, B. T., De Silva, G. M., & Lattanzio, J. 2015, MNRAS, 446, 3556, doi: 10.1093/mnras/stu2348
  * Marino (2013) Marino, A. F. 2013, Mem. Soc. Astron. Italiana, 84, 29, doi: 10.48550/arXiv.1302.2306
  * Marino et al. (2011) Marino, A. F., Villanova, S., Milone, A. P., et al. 2011, ApJ, 730, L16, doi: 10.1088/2041-8205/730/2/L16
  * Marino et al. (2012) Marino, A. F., Milone, A. P., Piotto, G., et al. 2012, The Astrophysical Journal, 746, 14, doi: 10.1088/0004-637X/746/1/14
  * Marino et al. (2016) Marino, A. F., Milone, A. P., Casagrande, L., et al. 2016, MNRAS, 459, 610, doi: 10.1093/mnras/stw611
  * Marley et al. (2002) Marley, M. S., Seager, S., Saumon, D., et al. 2002, ApJ, 568, 335, doi: 10.1086/338800
  * Martell & Grebel (2010) Martell, S. L., & Grebel, E. K. 2010, A&A, 519, A14, doi: 10.1051/0004-6361/201014135
  * Martell & Smith (2009) Martell, S. L., & Smith, G. H. 2009, PASP, 121, 577, doi: 10.1086/599979
  * Martell et al. (2011) Martell, S. L., Smolinski, J. P., Beers, T. C., & Grebel, E. K. 2011, A&A, 534, A136, doi: 10.1051/0004-6361/201117644
  * Mashonkina (2010) Mashonkina, L. 2010, in EAS Publications Series, Vol. 43, EAS Publications Series, ed. R. Monier, B. Smalley, G. Wahlgren, & P. Stee, 189–197, doi: 10.1051/eas/1043014
  * Messenger & Lattanzio (2002) Messenger, B. B., & Lattanzio, J. C. 2002, MNRAS, 331, 684, doi: 10.1046/j.1365-8711.2002.05234.x
  * Messier (1781) Messier, C. 1781, Catalogue des Nébuleuses et des Amas d’Étoiles (Catalog of Nebulae and Star Clusters), Connoissance des Temps ou des Mouvements Célestes, for 1784, p. 227-267
  * Mészáros & Allende Prieto (2013) Mészáros, S., & Allende Prieto, C. 2013, MNRAS, 430, 3285, doi: 10.1093/mnras/stt130
  * Mészáros et al. (2012) Mészáros, S., Allende Prieto, C., Edvardsson, B., et al. 2012, AJ, 144, 120, doi: 10.1088/0004-6256/144/4/120
  * Mihalas (1978) Mihalas, D. 1978, Stellar atmospheres (W. H. Freeman)
  * Milone (2023) Milone, A. P. 2023, arXiv e-prints, arXiv:2304.07770, doi: 10.48550/arXiv.2304.07770
  * Milone & Marino (2022) Milone, A. P., & Marino, A. F. 2022, Universe, 8, 359, doi: 10.3390/universe8070359
  * Milone et al. (2013) Milone, A. P., Marino, A. F., Piotto, G., et al. 2013, ApJ, 767, 120, doi: 10.1088/0004-637X/767/2/120
  * Milone et al. (2019) Milone, A. P., Marino, A. F., Bedin, L. R., et al. 2019, MNRAS, 484, 4046, doi: 10.1093/mnras/stz277
  * Moehler (2001) Moehler, S. 2001, PASP, 113, 1162, doi: 10.1086/323297
  * Nardiello et al. (2023) Nardiello, D., Griggio, M., & Bedin, L. R. 2023, MNRAS, 521, L39, doi: 10.1093/mnrasl/slad021
  * Nendwich et al. (2004) Nendwich, J., Heiter, U., Kupka, F., Nesvacil, N., & Weiss, W. W. 2004, Communications in Asteroseismology, 144, 43, doi: 10.48550/arXiv.astro-ph/0406381
  * Newberg (2016) Newberg, H. J. 2016, in Astrophysics and Space Science Library, Vol. 420, Tidal Streams in the Local Group and Beyond, ed. H. J. Newberg & J. L. Carlin, 1, doi: 10.1007/978-3-319-19336-6_1
  * Nitschai et al. (2023) Nitschai, M. S., Neumayer, N., Clontz, C., et al. 2023, arXiv e-prints, arXiv:2309.02503, doi: 10.48550/arXiv.2309.02503
  * Norris (1981) Norris, J. 1981, ApJ, 248, 177, doi: 10.1086/159141
  * Norris & Bessell (1975) Norris, J., & Bessell, M. S. 1975, ApJ, 201, L75, doi: 10.1086/181946
  * Norris et al. (1981) Norris, J., Cottrell, P. L., Freeman, K. C., & Da Costa, G. S. 1981, ApJ, 244, 205, doi: 10.1086/158698
  * Norris & Freeman (1979) Norris, J., & Freeman, K. C. 1979, ApJ, 230, L179, doi: 10.1086/182988
  * Oort (1926) Oort, J. H. 1926, PhD thesis, -
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Percival et al. (2002) Percival, S. M., Salaris, M., van Wyk, F., & Kilkenny, D. 2002, ApJ, 573, 174, doi: 10.1086/340593
  * Peterson (1980) Peterson, R. C. 1980, ApJ, 237, L87, doi: 10.1086/157950
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38, doi: 10.1051/0004-6361/201937381
  * Piotto et al. (2007) Piotto, G., Bedin, L. R., Anderson, J., et al. 2007, ApJ, 661, L53, doi: 10.1086/518503
  * Piotto et al. (2012) Piotto, G., Milone, A. P., Anderson, J., et al. 2012, ApJ, 760, 39, doi: 10.1088/0004-637X/760/1/39
  * Plez (2011) Plez, B. 2011, in Journal of Physics Conference Series, Vol. 328, Journal of Physics Conference Series, 012005, doi: 10.1088/1742-6596/328/1/012005
  * Popper (1947) Popper, D. M. 1947, ApJ, 105, 204, doi: 10.1086/144893
  * Qin et al. (2023) Qin, S., Zhong, J., Tang, T., & Chen, L. 2023, ApJS, 265, 12, doi: 10.3847/1538-4365/acadd6
  * Ramírez et al. (2012) Ramírez, I., Meléndez, J., & Chanamé, J. 2012, ApJ, 757, 164, doi: 10.1088/0004-637X/757/2/164
  * Recio-Blanco (2018) Recio-Blanco, A. 2018, A&A, 620, A194, doi: 10.1051/0004-6361/201833179
  * Rennó et al. (2020) Rennó, C., Barbuy, B., Moura, T. C., & Trevisan, M. 2020, MNRAS, 498, 5834, doi: 10.1093/mnras/staa2697
  * Ribas et al. (2015) Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52, doi: 10.1051/0004-6361/201424846
  * Rood & Crocker (1985) Rood, R. T., & Crocker, D. A. 1985, in Horizontal-Branch and UV-Bright Stars, ed. A. G. D. Philip, 99–105
  * Salaris & Cassisi (2005) Salaris, M., & Cassisi, S. 2005, Evolution of Stars and Stellar Populations (John Wiley and Sons)
  * Salaris & Cassisi (2014) —. 2014, A&A, 566, A109, doi: 10.1051/0004-6361/201423722
  * Salaris & Cassisi (2015) —. 2015, A&A, 577, A60, doi: 10.1051/0004-6361/201525812
  * Salaris et al. (2020) Salaris, M., Usher, C., Martocchia, S., et al. 2020, MNRAS, 492, 3459, doi: 10.1093/mnras/staa089
  * Sampedro et al. (2017) Sampedro, L., Dias, W. S., Alfaro, E. J., Monteiro, H., & Molino, A. 2017, MNRAS, 470, 3937, doi: 10.1093/mnras/stx1485
  * Sandage & Wildey (1967) Sandage, A., & Wildey, R. 1967, ApJ, 150, 469, doi: 10.1086/149350
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713, doi: 10.1086/192204
  * Sbordone & Bonifacio (2005) Sbordone, L., & Bonifacio, P. 2005, Atlas Cookbook.  http://atmos.obspm.fr/index.php/documentation/7
  * Smith (1987) Smith, G. H. 1987, PASP, 99, 67, doi: 10.1086/131958
  * Smith & Norris (1982a) Smith, G. H., & Norris, J. 1982a, ApJ, 254, 149, doi: 10.1086/159718
  * Smith & Norris (1982b) —. 1982b, ApJ, 254, 594, doi: 10.1086/159770
  * Smith & Norris (1983) —. 1983, ApJ, 264, 215, doi: 10.1086/160588
  * Sneden et al. (1992) Sneden, C., Kraft, R. P., Prosser, C. F., & Langer, G. E. 1992, AJ, 104, 2121, doi: 10.1086/116388
  * Sollima & Baumgardt (2017) Sollima, A., & Baumgardt, H. 2017, MNRAS, 471, 3668, doi: 10.1093/mnras/stx1856
  * Spada et al. (2017) Spada, F., Demarque, P., Kim, Y. C., Boyajian, T. S., & Brewer, J. M. 2017, ApJ, 838, 161, doi: 10.3847/1538-4357/aa661d
  * Strader et al. (2004) Strader, J., Brodie, J. P., & Forbes, D. A. 2004, AJ, 127, 3431, doi: 10.1086/420995
  * Sweigart & Mengel (1979) Sweigart, A. V., & Mengel, J. G. 1979, ApJ, 229, 624, doi: 10.1086/156996
  * ten Bruggencate (1927) ten Bruggencate, P. 1927, Sternhaufen: Ihr Bau, Ihre Stellung zum Sternsystem und Ihre Bedeutung für die Kosmogonie (Springer), 127–135
  * Thompson et al. (2020) Thompson, I. B., Udalski, A., Dotter, A., et al. 2020, MNRAS, 492, 4254, doi: 10.1093/mnras/staa032
  * Thygesen et al. (2014) Thygesen, A. O., Sbordone, L., Andrievsky, S., et al. 2014, A&A, 572, A108, doi: 10.1051/0004-6361/201424533
  * Tsuji & Nakajima (2003) Tsuji, T., & Nakajima, T. 2003, ApJ, 585, L151, doi: 10.1086/374388
  * Valcin et al. (2020) Valcin, D., Bernal, J. L., Jimenez, R., Verde, L., & Wandelt, B. D. 2020, J. Cosmology Astropart. Phys, 2020, 002, doi: 10.1088/1475-7516/2020/12/002
  * Valcin et al. (2021) Valcin, D., Jimenez, R., Verde, L., Bernal, J. L., & Wandelt, B. D. 2021, J. Cosmology Astropart. Phys, 2021, 017, doi: 10.1088/1475-7516/2021/08/017
  * van den Bergh (1967) van den Bergh, S. 1967, AJ, 72, 70, doi: 10.1086/110203
  * VandenBerg (2000) VandenBerg, D. A. 2000, ApJS, 129, 315, doi: 10.1086/313404
  * VandenBerg et al. (2013) VandenBerg, D. A., Brogaard, K., Leaman, R., & Casagrande, L. 2013, ApJ, 775, 134, doi: 10.1088/0004-637X/775/2/134
  * Wallerstein (1969) Wallerstein, G. 1969, ApJ, 158, 607, doi: 10.1086/150222
  * Wallerstein & Greenstein (1964) Wallerstein, G., & Greenstein, J. L. 1964, ApJ, 139, 1163, doi: 10.1086/147858
  * Wenger & Champion (1998) Wenger, C., & Champion, J. P. 1998, J. Quant. Spec. Radiat. Transf., 59, 471, doi: 10.1016/S0022-4073(97)00106-4
  * Wing & Stock (1973) Wing, R. F., & Stock, J. 1973, ApJ, 186, 979, doi: 10.1086/152561
  * Winter & Clarke (2023) Winter, A. J., & Clarke, C. J. 2023, MNRAS, 521, 1646, doi: 10.1093/mnras/stad312
  * Yong et al. (2008) Yong, D., Meléndez, J., Cunha, K., et al. 2008, ApJ, 689, 1020, doi: 10.1086/592229
  * Zhang et al. (2008) Zhang, H. W., Gehren, T., & Zhao, G. 2008, A&A, 481, 489, doi: 10.1051/0004-6361:20078910
  * Zhou et al. (2022) Zhou, T., Whisnant, A., Gerasimov, R., & Burgasser, A. J. 2022, Research Notes of the American Astronomical Society, 6, 212, doi: 10.3847/2515-5172/ac9ab6
  * Ziliotto et al. (2023) Ziliotto, T., Milone, A., Marino, A. F., et al. 2023, ApJ, 953, 62, doi: 10.3847/1538-4357/acde76
  * Zinn (1973) Zinn, R. 1973, A&A, 25, 409