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

Stellar evolution is simulated using the Modules for Experiments in Stellar Astrophysics (MESA) 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 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 ε_nuc > 10³ erg g⁻¹ s⁻¹ begins to fall outside the central 10% of the stellar mass (ε_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|{\rm{\Delta }}{\mathrm{log}}_{10}L\right|\lt 0.12\,\mathrm{dex}$ and $\left|{\rm{\Delta }}{T}_{\mathrm{eff}}\right|\lt 120\,{\rm{K}}$, respectively, at the checkpoint ages between 10 and 13.5 Gyr in steps of 0.5 Gyr (the typical average differences at 13.5 Gyr are $\langle \left|{\rm{\Delta }}{\mathrm{log}}_{10}L\right|\rangle \approx 0.07\,\mathrm{dex}$ and $\langle \left|{\rm{\Delta }}{T}_{\mathrm{eff}}\right|\rangle \approx 70\,{\rm{K}}$). The lowest mass in the grid (∼0.04–0.05 M_⊙) is chosen to attain an evolved T_eff of ≲700 K at 13.5 Gyr, while the highest mass (∼0.9 M_⊙) is set by requiring that the star take at least 10 Gyr to reach the tip of the subgiant branch.

PMS stars at zero age are initialized with uniform chemical abundances and a central temperature of T_c = 5e5 K. 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 >0.3 M_⊙ 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 calculations, ⁸ and it 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 α_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 website ⁹ and Zenodo doi:10.5281/zenodo.10016008.

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 owing to the complex relationship between the pressure–temperature profile and the outgoing energy flux, caused by low-temperature phenomena such as nongray molecular opacity, cloud formation, and convective instabilities. This issue is particularly prominent at T_eff ≲ 4500 K, 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\,\leqslant {\mathrm{log}}_{10}(g)\leqslant 6$ for 2500 K ≤ T_eff ≤ 7500 K and $4\leqslant {\mathrm{log}}_{10}(g)\leqslant 6$ for 500 K ≤ T_eff ≤ 2400 K, with steps of 100 K in T_eff and 0.5 in ${\mathrm{log}}_{10}(g)$. Models with ${\mathrm{log}}_{10}(g)\lt 4\ \mathrm{and}\ {T}_{\mathrm{eff}}\lt 2500\,{\rm{K}}$ were not required, since such low gravities are only encountered during the early evolution (<5 Myr) and the subgiant phase, both of which are characterized by higher effective temperatures (T_eff > 2500 K). For some of the calculated isochrones, a small number of model atmospheres could not be computed owing 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_eff–${\mathrm{log}}_{10}(g)$ tracks of selected evolutionary models.

Zoom In Zoom Out Reset image size

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 (τ) 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 nonideal gas behavior may be significant (Burke et al. 2004), and is not accounted for in our high-T_eff model atmospheres. We therefore employed two distinct atmosphere–interior coupling regimes with boundary conditions at τ = 100 for M < 0.5 M_⊙, and at T(τ) = T_eff (photosphere) for M > 0.6 M_⊙, where T(τ) is the temperature profile of the star and M is the initial stellar mass. At intermediate stellar masses, 0.5 M_⊙ ≤ M ≤ 0.6 M_⊙, a linear ramp between the two regimes was applied. The range of transition was chosen to approximately coincide with the maximum value of dT(τ = 100)/dM for main-sequence stars. The implications of this choice for the synthetic photometry are discussed in Section 4.6.

When the range of pretabulated 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_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 website ¹⁰ and Zenodo. ¹¹

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_eff ≤ 5000 K and the ATLAS code (Kurucz 1970), version 9 (Kurucz 2005, 2014), at T_eff > 5000 K. The use of ATLAS for "warm" atmospheres is advantageous owing to its superior computational efficiency (see Section 4.5), attained by virtue of opacity sampling from pretabulated 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 9 nm ≲ λ ≲ 160 μm at 57 temperatures between ≈2e3 K and ≈2e5 K using the satellite package DFSYNHTE (Castelli & Kurucz 2003; Castelli 2005). ATLAS atmospheres are stratified into 72 plane-parallel layers, spanning ∼10⁻⁷ < τ < 10¹⁰⁰ with logarithmic spacing. The synthetic spectra for each model atmosphere were calculated using the SYNTHE code (Kurucz & Avrett 1981) within 0.1 μm < λ < 5 μm, at the constant resolution ¹² of λ/Δλ ∼ 6 × 10⁵. All three codes are packaged in the open-source Python dispatcher 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_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₂O lines from Barber et al. 2006 and CH₄ 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 nonideal 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_eff ≲ 3000 K compared to gas-only models (Lunine et al. 1986; Allard et al. 2001). Our setup of PHOENIX also includes the Allard & Homeier treatment of gravitational settling for condensate species (Helling et al. 2008; Allard et al. 2012) that gradually removes clouds from the spectrum-forming region of the atmosphere at T_eff ≪ 2000 K (Tsuji & Nakajima 2003). For computational efficiency, we only consider gravitational settling at T_eff < 2500 K 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 (τ_λ , where λ = 1.275 μm) instead of the Rosseland mean. ¹⁴ The bottom of the atmosphere is set to τ_λ = 10³ for the models without gravitational settling but reduced to τ_λ = 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 10⁻³ dyn cm⁻². 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 τ = 100 layer is always found within the outer 0.1% of the photospheric radius, calculated at T(τ) = T_eff). Therefore, the effect of spherical geometry is expected to be subdominant, allowing the use of approximate atmospheric radii from the 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 1 Å < λ < 1 mm with a median resolution of λ/Δλ ∼ 18,250 between 0.5 and 3 μm. A few examples of the synthetic spectra for the chemical abundances of the ridgeline isochrone and ${\mathrm{log}}_{10}(g)=5.0$ are provided in Figure 2. The atmosphere models calculated in this study are made available online. ¹⁵

Zoom In Zoom Out Reset image size

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], which is reduced by 3 dex in all PHOENIX models compared to the PMS abundance. This assumption is justified by the fact that low-mass stars (≲0.6 M_⊙) 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).

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_eff > 6500 K 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_eff PHOENIX atmospheres is more involved. In models with subdominant cloud settling (T_eff ≳ 2000 K), 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 ∼3% and 0.5%, respectively. 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_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_eff and ${\mathrm{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 owing to the rapid growth of the oscillation amplitude in the temperature structure between iterations. For the rest of the 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 ∼100 iterations. The models with limited convergence are shown in Figure 1 with red markers for the ridgeline isochrone. At very low temperatures (T_eff ≲ 1200 K) and high gravities (${\mathrm{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_eff and ${\mathrm{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 unparallelized 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 (≲0.05 M_⊙), where the tolerances were relaxed to their standard values to allow evolution over the discontinuities in ${dT}(\tau =100)/d{\mathrm{log}}_{10}(g)$, caused by the sparse gravity sampling in the boundary condition tables.

All model atmospheres used in this work were calculated on the Bridges2 supercomputer (Brown et al. 2021), operated by the Pittsburgh Supercomputing Center. Our PHOENIX models were calculated to ∼100 iterations, requiring between ∼0.15 processor hours per iteration at the highest temperatures to ∼1.5 processor hours at the lowest. The high memory demand of PHOENIX (≳4 GB per model) necessitated requesting at least two (and occasionally three) 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 ≈0.005 processor hours per model for all iterations.

The computational demand of a complete atmosphere grid for each chemical composition is approximately 5 × 10⁴ 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.

The process of calculating theoretical isochrones described earlier in this section required a choice of three "stitching points," where distinct modeling setups were blended together: the effective temperature where ATLAS models are replaced with PHOENIX models (5000 K), the effective temperature where the gravitational settling in PHOENIX models was disabled (2500 K), and the initial stellar mass, where the τ = 100 boundary condition tables were replaced with T(τ) = T_eff (smooth ramp between 0.5 and 0.6 M_⊙). Here we briefly review the implications of those choices for synthetic photometry.

Our choice of the ATLAS/PHOENIX transition temperature (5000 K) is higher than that used in our previous work (4000 K; Gerasimov et al. 2022a). For comparison, the published grid of ATLAS models with the updated molecular opacities (Castelli & Kurucz 2003) reaches 3500 K, 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 4000 K ≤ T_eff ≤ 5000 K, as well as an additional grid of evolutionary models with the updated boundary conditions. We found the synthetic photometry at 11.5 Gyr of the regular ridgeline isochrone and its altered version to be indistinguishable at T_eff = 5000 K but rapidly diverging at lower temperatures to attain the difference of ≈0.025 mag in F606W − F814W and F150W2 − F322W2 at T_eff = 4000 K. The difference originates almost entirely from the synthetic spectra, with only a minor contribution from the boundary conditions.

Our choice of the minimum T_eff without gravitational settling (2500 K) is, on the other hand, lower than that employed in our previous study of ω Centauri (3000 K; 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 2500 K ≤ T_eff ≤ 3000 K and enabled settling of condensates, as well as a test set of corresponding evolutionary models. We found that at T_eff = 2500 K and at the age of 11.5 Gyr the effect of gravitational settling on the same colors as before is ≈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_eff > 2500 K, 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 3000 K, 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 0.5 M_⊙ ≤ M ≤ 0.6 M_⊙. The effect on synthetic photometry at 11.5 Gyr is most prominent in F606W − F814W when the range of T(τ) = T_eff boundary condition tables is extended down to 0.5 M_⊙, resulting in the discontinuity of ≲0.01 mag between the two regimes. The discontinuity is less prominent when the τ = 100 tables are extended to 0.6 M_⊙ 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 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_eff > 2500 K.

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_eff and ${\mathrm{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 13.2418 ± 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_eff spectra on surface parameters and the highly nonlinear 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 (Gustafsson et al. 2008; Mészáros & Allende Prieto 2013). 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 nonconvergent 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 600 K ≤ T_eff < 2400 K in 50 K steps, with ${\mathrm{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_eff–${\mathrm{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_eff < 2500 K from higher-gravity model atmospheres that have better convergence (e.g., from Figure 1, to calculate the bolometric corrections at the T_eff = 2000 K point on the 11.5 Gyr isochrone, one would require both ${\mathrm{log}}_{10}(g)\,=5.0$ and ${\mathrm{log}}_{10}(g)=5.5$ model atmospheres for the interpolation method, but only one ${\mathrm{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_eff–${\mathrm{log}}_{10}(g)$ relationship.

The linear interpolation and hammering methods of calculating synthetic photometry at T_eff ≤ 2350 K are compared in Figure 3. The figure demonstrates that at T_eff ≲ 1000 K (the brown dwarf regime) linear interpolation may be off by as much as 2 mag in the optical regime and ≳0.5 mag in the near-infrared. It is also apparent that the two methods converge at T_eff ≳ 2000 K, 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_eff > 2350 K.

Zoom In Zoom Out Reset image size

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_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.

The high effective temperatures of the main-sequence turnoff and the subgiant branch in 47 Tucanae (T_eff ≳ 5000 K) 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, chap. 5.5). These high-mass stars are also sensitive to the adopted energy transport parameters (in particular, the mixing length, α_MLT) owing 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 turnoff and the subgiant branch before analyzing the lower main sequence. We used this part of the CMD to assess our choice of α_MLT = 1.82; to fix the cluster age and the optical reddening, E(B − V); and to evaluate the accuracy of the helium mass fraction and metallicity in the nominal chemical composition (Section 2, Y = 0.25, [Fe/H] = − 0.75).

The top panels of Figure 4 show the turnoff point and the subgiant branch of the ridgeline isochrone, evaluated at four distinct ages and overplotted on the CMD 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) = 0.04 mag, in agreement with Percival et al. (2002). In the figure, the 11.5 Gyr isochrone is the only one that matches both the color of the turnoff point and the luminosity of the subgiant branch within one standard deviation. Furthermore, the presented fit firmly fixes the reddening at 0.04 mag, as lower values would result in the isochrone being "too blue" at the turnoff 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; see also Rennó et al. 2020, for a compilation of recent age measurements).

Zoom In Zoom Out Reset image size

To explore the effect of the helium mass fraction and the overall metallicity, we calculated two additional isochrones with Y = 0.28 and [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 turnoff 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 turnoff 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 bottom panels of the figure demonstrate how our choice of the mixing length compares to two alternative values: α_MLT = 1.6 and α_MLT = 2.0. The isochrones are shown at 11.5 Gyr. Around the turnoff point, the effect is practically indistinguishable from that of the cluster age, emphasizing the limitations of using the color of the turnoff point as an age diagnostic for globular clusters. Our earlier claim that the 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 12.5 Gyr for α_MLT = 2.0, or 10.5 Gyr for α_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 α_MLT = 1.82 for all isochrones and the corresponding best-fit age of 11.5 Gyr.

Below the main-sequence inflection point (∼0.5 M_⊙; 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 endpoint of the evolutionary track to different T_eff and ${\mathrm{log}}_{10}(g)$. The importance of the latter effect approximately correlates with $\delta \kappa =\left|d\kappa /d[{\rm{X}}/\mathrm{Fe}]\right|$ at T(τ) = T_eff, where [X/Fe] is the abundance of the element of interest and κ is the Rosseland mean opacity. For the range 0.1 M_⊙ < M < 0.5 M_⊙, this diagnostic is by far the largest for [O/Fe] and [Ti/Fe] (∼ 10⁻⁷ cm⁻¹ dex⁻¹ at 0.2 M_⊙) owing to the prominent TiO and H₂O absorption bands (see Figure 2). For comparison, the next two most important elements, [C/Fe] and [Al/Fe], have δκ ∼ 2 × 10⁻⁸ cm⁻¹dex⁻¹ and ∼10⁻⁸ cm⁻¹ dex⁻¹, 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_eff and ${\mathrm{log}}_{10}(g)$ are only weakly sensitive to the oxygen abundance, since H₂O 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 δ κ, and since both are considered to be α-elements, even solar-scaled boundary conditions with adjustable α-enhancement would not be adequate.

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 data sets only extend to the end of the main sequence (Section 3), we terminate all intermediate isochrones at T_eff = 3000 K to avoid unnecessary calculations of model atmospheres.

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 five initial stellar masses between 0.13 and 0.7 M_⊙ 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_eff and ${\mathrm{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 +0.5 and −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 δ κ for each element to estimate our confidence in the derived correlation between colors and abundances.

Zoom In Zoom Out Reset image size

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 δ κ 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_eff ≳ 700 K to reach the substellar regime, as required for our further analysis. For the final ridgeline isochrone, an additional abundance variation grid was computed for seven initial masses below the end of the main sequence at 0.05 M_⊙ ≤ M ≤ 0.08 M_⊙. 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, and 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 ∼19 mag). A similar red excess in the optical CMD is also seen at the corresponding evolutionary phase (F814W ∼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.

Zoom In Zoom Out Reset image size

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 ([O/Fe] = 0.48 dex for the blue tail isochrone, [O/Fe] = 0.33 dex for ridgeline, and [O/Fe] = 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–[O/Fe] space.

To obtain the desired transformation, we first generated a 1000 × 300 regular grid of synthetic stars at 1000 evenly spaced initial masses from the minimum (≈0.045 M_⊙) to the maximum (≈0.875 M_⊙) value within the range of all three isochrones and 300 evenly spaced values of [O/Fe] from −0.57 to 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 nonunique. We carry out this analysis below the knee (0.1 M_⊙ ≤ M ≤ 0.5 M_⊙), where the isochrones are well separated, and the transformation is accurate to 10⁻³ M_⊙ in mass and 10⁻² dex in oxygen abundance, as estimated by applying the transformation to 10⁴ additional synthetic stars with randomly chosen masses and oxygen abundances. We further restrict our analysis to stars that satisfy −0.27 dex < [O/Fe] < 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.

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

Zoom In Zoom Out Reset image size

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⁵ synthetic equivalents of the spectroscopic data set 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, and 5th, 25th, 75th, and 95th percentiles. The averages and spreads in those statistics across all 10⁵ Monte Carlo trials are summarized in Table 2, alongside the equivalent values calculated for the real spectroscopic data set. 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 ∼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 ≈0.03 dex at the oxygen-poor tail of the distribution (inferred from the 5th percentile) to ≈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 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 bottom 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 ± 0.03 and −0.43 ± 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. ∼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.

The mass functions of globular clusters retain a footprint of their dynamical evolution (Baumgardt & Sollima 2017; Sollima & Baumgardt 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 subsection, 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 (ξ) is of the form of a broken power law (Kroupa 2001)

$$\begin{eqnarray}\xi (M)\propto \left\{\begin{array}{ll}{M}^{-\beta }, & \mathrm{if}\ M\gt {M}_{\mathrm{bp}}\\ {M}^{-\gamma }, & \mathrm{if}\ M\leqslant {M}_{\mathrm{bp}}.\end{array}\right.\end{eqnarray} \tag{ 2 }$$

Here β and γ are the power-law indices in the high- and low-mass regimes, respectively, separated by the break-point stellar mass M_bp. We keep all three parameters as free variables. The magnitude function, ϕ(m), where m is the apparent magnitude of the star (in our case, F160W), is related to the mass function as

$$\begin{eqnarray}&&\phi (m)=-\xi \left(M(m)\right)\displaystyle \frac{{dM}(m)}{{dm}}.\end{eqnarray} \tag{ 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:

$$\begin{eqnarray}&&\phi \left(\displaystyle \frac{{m}_{i}+{m}_{i+1}}{2}\right)=-\xi \left(\displaystyle \frac{{M}_{i}+{M}_{i+1}}{2}\right)\displaystyle \frac{{M}_{i+1}-{M}_{i}}{{m}_{i+1}-{m}_{i}}.\end{eqnarray} \tag{ 4 }$$

Here the magnitude function is evaluated at the midpoints 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 0.15 M_⊙ < M < 0.4 M_⊙ and −0.27 dex < [O/Fe] < 0.63 dex (the same restrictions as in the bottom 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 owing to contamination by SMC members and potential incompleteness of the photometric sample.

To determine the parameters of the mass function, β, γ, and M_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 β = 1.75 ± 0.15, γ = 0.08 ± 0.04, and M_bp = 0.47 ± 0.02 M_⊙. 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.

Zoom In Zoom Out Reset image size

In this subsection, 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 versus F150W2 − F322W2, since brown dwarfs are most likely to be detected in wide bands at long wavelengths owing 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 recalculated the synthetic magnitude function of 47 Tucanae in 30 evenly spaced magnitude bins within the range of all three theoretical isochrones, using the best-fit mass function parameters. The calculations were carried out for a variety of ages between 0.1 and 13.5 Gyr, including 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 ∼ 26, where the brown dwarf density per magnitude is reduced by nearly an order of magnitude (compared to the maximum density within the modeled range that is predicted to occur around m = 27.5). While the lower main sequence is virtually unaffected by the cluster age, both the turnoff 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 11.5 Gyr is shown in Figure 10. Lines of constant mass (isomass lines) for selected masses and −0.27 dex < [O/Fe] < 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.

Zoom In Zoom Out Reset image size

The predicted mass–magnitude relationship for the ridgeline isochrone is shown in an inset in 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 11.5 Gyr. In this case, the HBL was computed as 0.074 M_⊙.