© The Authors 2025

1. Introduction

Globular clusters (GCs), comprising the oldest and largest coeval stellar populations, serve as crucial markers for the chemical evolution of the Galaxy and provide valuable insights into the nature of ancient metal-poor stars. The long-standing belief that all stars within a cluster have the same chemical composition was challenged by spectroscopic studies, which revealed the presence of chemical inhomogeneities in some clusters (see e.g. Gratton et al. 2012, and references therein). Two decades of space-based photometric studies have further revolutionised our understanding of GCs, uncovering intricate colour-magnitude diagrams (CMDs) that reveal the presence of multiple stellar populations (mPOPs) with varying chemical compositions. Despite these discoveries, the process driving this chemical diversity remains unclear, partly because most research on GCs has focused on the brighter stars, such as those on the upper main sequence (MS) or in the giant phase (see Bastian & Lardo 2018; Gratton et al. 2019, for a recent review).

In this paper, we present JWST images of NGC 5139 (hereafter referred to as ω Cen), the relatively close (we assumed the distance of 5.24 ± 0.11 kpc from Soltis et al. 2021, throughout the paper) and most massive (∼4 × 10⁶ M_⊙; Giersz & Heggie 2003; D’Souza & Rix 2013) GC in the Milky Way. These observations aim to investigate the faintest objects within the cluster using near-infrared (NIR) data and explore the mPOPs in the faintest regions of its CMD.

ω Cen is well known for its complex mPOP system, making it one of the most enigmatic stellar systems in the Galaxy. It hosts at least two main stellar components: the blue main sequence (bMS) and the red main sequence (rMS) (e.g. Bedin et al. 2004), with the bMS more centrally concentrated than the rMS (Bellini et al. 2009; Sollima et al. 2007a). Except for the cluster’s central region, the rMS is the dominant population. This initially appeared inconsistent with spectroscopic studies, which indicated that lower-metallicity stars are the most abundant in the cluster (e.g. Norris & Da Costa 1995). To resolve this discrepancy, Bedin et al. (2004) proposed that the bMS might be significantly more enhanced in helium than the rMS. This hypothesis was supported by Norris (2004), who estimated a helium mass fraction difference of ΔY ∼ 0.12 between the two sequences. Further confirmation came from Piotto et al. (2005), who demonstrated that bMS stars are indeed more metal-rich than their rMS counterparts by ∼0.3 dex, while King et al. (2012) derived a helium mass fraction of Y = 0.39 ± 0.02 for the bMS. However, a more recent spectroscopic analysis by Latour et al. (2021) based on MUSE data suggests a smaller helium difference between the bMS and rMS (Y ≲ 0.1).

Another remarkable feature of ω Cen is the unusually broad metallicity spread observed within each sequence, significantly larger than what is typically found in other GCs. The metallicity ranges from [Fe/H] ≈ −2.2 to −0.5, with [X/Fe] dispersions reaching 0.5 to more than 1.0 dex for several elements (e.g. Norris & Da Costa 1995; Suntzeff & Kraft 1996; Smith et al. 2000; Johnson et al. 2008, 2009; Bellini et al. 2017a). This suggests that ω Cen might be the nucleus of a dwarf galaxy absorbed by the Milky Way, or the outcome of the merger of two or more clusters (Norris et al. 1997; Jurcsik 1998; Bekki & Freeman 2003; Pancino et al. 2000; Bekki & Norris 2006; Ibata et al. 2019; van de Ven et al. 2006). Recent studies have further revealed that both the bMS and rMS are subdivided into multiple subpopulations, with up to 15 distinct stellar populations identified within the cluster so far (Bellini et al. 2017a; Scalco et al. 2024a).

The high-resolution deep NIR images provided by JWST have proven to be very effective in the exploration of mPOPs in GCs (47 Tuc; Milone et al. 2023; Marino et al. 2024, M92; Nardiello et al. 2022, NGC 6440; Cadelano et al. 2023 and NGC 6397; Scalco et al. 2024b). Here, we present the first¹ analysis of the mPOPs in ω Cen using JWST photometry, with a focus on their luminosity function (LF) and mass function (MF) along the MS.

The paper is organised as follows. Section 2 describes the data and reduction methods, while Section 3 focuses on the evaluation of proper motions (PMs). Section 4 presents the CMD of ω Cen based on JWST photometry, followed by Section 5, which details the artificial star (AS) tests. Section 6 analyses the LFs and MFs of ω Cen and its stellar populations using JWST data. Finally, a brief summary is provided in Section 7.

2. Observations and data reduction

The data analysed in this study were obtained from two programmes:

• Our proprietary images from JWST GO-5110 programme (Bedin et al. 2024a). The JWST data were obtained with the Near Infrared Camera (NIRCam, Rieke et al. 2023) simultaneously with the Short Wavelength (SW) and Long Wavelength (LW) channels on August 7–8, 2024 (epoch ∼2024.6). Ultra-wide filters were used for both NIRCam channels: F150W2 for the SW channel and F322W2 for the LW channel. The six-point FULLBOX primary pattern was employed with a THREE-POINT-LARGE-WITH-NIRISS sub-pixel dither pattern. At each of the resulting 18 pointings, a single image was captured in both channels using the MEDIUM8 readout pattern (nine groups), resulting in an effective exposure time of 944.836 s per image.

• Archival material from our HST GO-14118+14662 (Bedin et al. 2016a,b) multi-epoch programme. This dataset includes images taken with the Wide Field Channel (WFC) of the Advanced Camera for Surveys (ACS) using the F606W and F814W filters, spanning a period from ∼2015.6 to ∼2018.5. These data were reduced and presented in a recent paper by Scalco et al. (2024c), where a comprehensive description of the dataset can be found.

Figure 1 shows the positions of our datasets within the field of view (FOV), overlaid on a Digital Sky Survey² (DSS) image of ω Cen. The overlapping region between the two datasets, highlighted in magenta, covers a radial range from ∼2.4 r_h to ∼3.2 r_h.

Fig. 1. NIRCam JWST field from the GO-5110 programme (red) and the ACS/WFC HST field from the GO-14118+14662 programme (blue) overlaid on a DSS image of ω Cen. The overlap region between the two datasets is highlighted in magenta. Units are in arcminutes from the cluster’s centre. The dashed white and green circles represent the core radius (rc = 2′⋅37) and the half-light radius (rh = 5′⋅00), respectively, with additional circles marking 2 rh, 3 rh, and 4 rh. The values of rc and rh are from Trager et al. (1995), McLaughlin & van der Marel (2005) as reported in the Harris (1996, 2010) catalogue.

JWST images were processed using the software tools and methods detailed in Papers I, II, and III of the ‘Photometry and Astrometry with JWST’ series (Nardiello et al. 2022, 2023a; Griggio et al. 2023), which have been successfully applied in recent studies of 47 Tuc by Nardiello et al. (2023b), Scalco et al. (2025). We first processed the level-1b uncalibrated (_uncal) images using a development version of the JWST pipeline³ (Bushouse et al. 2023), running it through stages 1 and 2 to produce the level-2b (_cal) images. In the stage 1 pipeline, we utilised the default parameters with one exception: for the ramp fitting process, we employed the frame zero (the initial frame of each integration) to measure pixels that were saturated in the first group and extended up to the ramp. This approach effectively increased the dynamic range of our data by almost two magnitudes, extending into the nominally saturated intensity regime. The stage 2 pipeline was executed with all default parameters.

The reduction of the _cal images involved two main steps: first-pass and second-pass photometry. In the first-pass, starting from the empirical library point spread function (PSF) derived in Bedin et al. (2024b), we extracted PSFs, positions, and magnitudes of stars from each individual image. A geometric distortion correction was applied to the star positions using the solution by Griggio et al. (2023). The coordinates from each image were then transformed to a common reference frame, with bright cluster members from Gaia Data Release 3 (DR3; Gaia Collaboration 2016, 2023) serving as the reference, after transforming their positions to the epoch of the data collection.

For the second-pass photometry, we utilised a modified version of the KS2 code as described in Bellini et al. (2017b, 2018), Nardiello et al. (2018), Libralato et al. (2018, 2022), Scalco et al. (2021) and references therein. In this step, positions and fluxes were extracted using the PSFs and transformations obtained during the first-pass. KS2 processes all images simultaneously, making it well suited for detecting faint sources that may not be visible in individual frames. Along with fluxes and positions, KS2 generates several quality parameters, such as the PSF quality-of-fit (QFIT) parameter, the RADXS parameter (a metric for comparing the shape of a source to the PSF; see Bedin et al. 2008, 2009), and the local sky noise (rmsSKY). Detailed descriptions of these diagnostics can be found in Bellini et al. (2018), Scalco et al. (2021), Nardiello et al. (2018).

We applied a similar reduction process to the HST data, following the same first- and second-pass photometry approach. For a detailed explanation of the HST data reduction, we refer to Scalco et al. (2024c).

We calibrated our HST and JWST photometry to the VEGA-magnitude system following the procedures illustrated in Nardiello et al. (2023a, see also Bedin et al. 2005), while we anchored the astrometry to the International Celestial Reference System (ICRS) frame using Gaia DR3 data for sources in the observed fields.

We selected a sample of well-measured stars using the quality parameters provided by KS2 (QFIT, RADXS, and rmsSKY). In particular, we retained sources with QFIT > 0 and absolute RADXS values smaller than 0.05 in both filters. As for rmsSKY, we manually defined a fiducial threshold as a function of magnitude, following the trend of rmsSKY with magnitude, and excluded all sources lying above this threshold. Finally, we corrected our photometry for the effects of differential reddening and spatial zero-point variations following the procedure introduced in Sarajedini et al. (2007, see also Milone et al. 2012; Bellini et al. 2017c). Filter-specific extinction coefficients were computed assuming A_V = 3.1 × E(B − V), with a mean reddening E(B − V) = 0.12 (Harris 1996; Calamida et al. 2005). We adopted the extinction ratios for JWST/NIRCam filters from Wang & Chen (2019) (A_F150W2 = 0.15 × A_V, A_F322W2 = 0.04 × A_V). The measured differential reddening values δE(B − V) range from approximately −0.05 to +0.04 across the field. Figure 2 presents the m_F150W2 versus m_F150W2 − m_F322W2 CMD (panel a) and the m_F322W2 versus m_F150W2 − m_F322W2 CMD (panel b) for all stars in the JWST field. In both panels, black dots represent stars passing the photometric quality selections while grey dots represent stars not passing the selections. The black dots distinctly outline the cluster’s MS and low-MS, as well as its white dwarf (WD) cooling sequence (CS), located in the bottom-left of the CMDs. In both panels, the light red area indicates the onset of saturation, where photometry is obtained from the frame zero. The dark red region marks severe saturation, where even the frame zero is affected and the flux must be recovered from the unsaturated wings of the PSF. In this region, the photometry becomes unreliable, with the sequence displaying non-physical features. Both the photometry and the astrometry of stars in these regions must therefore be treated with particular caution. Stars within these saturated regions are not ideal for scientific analysis that require high photometric accuracy – such as the determination of the LFs and MFs of the individual stellar populations (see Section 6.2), which represents the core of this study – and will not be utilised for these purposes. However, for completeness and to provide a comprehensive view of the CMD morphology of ω Cen, we still include and make use of these stars in other parts of the paper. The reader should nevertheless keep in mind that any results involving stars in the saturated regime must be interpreted with appropriate caution.

Fig. 2. Colour-magnitude diagrams of ω Cen using JWST filters. (a) mF150W2 versus mF150W2 − mF322W2 CMD. (b) mF322W2 versus mF150W2 − mF322W2 CMD. In both panels, black dots represent stars passing the quality selections while grey dots represent stars not passing the selections. The shaded light red and dark red regions highlight areas of the CMDs affected by saturation. The light red region corresponds to saturated photometry, while the dark red region represents severe saturation, where the photometry is saturated even in the frame zero.

Figure 3 shows the same CMDs as in Fig. 2, but limited to stars located in the overlap region between the JWST and HST FOV. While the number of sources is noticeably reduced (by ∼30%), the two sequences along the MS remain visible. However, the presence of some contaminants highlights the need for a PM-based cluster membership selection, which will be addressed in the following section.

Fig. 3. Same as Fig. 2, but only for the sources located in the overlap region between the JWST and HST FOV.

3. Proper motions

Proper motions (PMs) were computed as displacements between the HST and JWST observations divided by the temporal baseline (∼7.5 years). For the HST epoch, we adopted the average position across the available exposures, given the much shorter time span relative to the HST–JWST baseline. The resulting PMs are shown in Fig. 4, which includes only stars that meet the photometric quality criteria described in the previous section and have measurable PMs. Panels (a) and (b) show the vector-point diagram (VPD) and the m_F322W2 versus m_F150W2 − m_F322W2 CMD, respectively. Since PMs are calculated relative to the cluster’s overall motion, the distribution of cluster members in the VPD is centred at (0,0). A second, sparser group of points is visible in the upper part of the VPD, representing background and foreground sources. Panel (c) shows the one-dimensional PM (μ_R, obtained by summing the PMs in the two directions in quadrature) plotted against m_F322W2. ω Cen members exhibit a narrow distribution in μ_R, mostly clustered below μ_R < 2 mas yr⁻¹, while field objects extend towards higher μ_R. We defined a conservative PM selection to separate cluster members from field objects, indicated by a red line. Panels (d) and (e) show the VPD and m_F322W2 versus m_F150W2 − m_F322W2 CMD for the stars that passed the PM selection, while panels (f) and (g) show the same diagrams for stars that did not pass the PM selection. In panels (b), (e), and (g), the shaded light red and dark red regions indicate areas of the CMDs affected by mild and severe saturation, respectively. The PMs of stars in these regions should be treated with caution. In what follows, our analysis will focus exclusively on the sample of stars shown in panels (d) and (e).

Fig. 4. Cluster membership selection based on PMs, for stars that passed the photometric quality selections and have measurable PMs. (a)–(b) VPD and mF322W2 versus mF150W2 − mF322W2 CMD, respectively. (c) mF322W2 magnitude versus the one-dimensional PM (μR). The red line separates cluster members from field stars. (d)–(e) VPD and mF322W2 versus mF150W2 − mF322W2 CMD for stars that passed the PM selection. (f)–(g) same diagrams for stars that did not pass the PM selection. In panels (c) through (g), sources that passed the PM selection are represented by blue dots, while those that did not pass are depicted as orange crosses. In panels (b), (e) and (g), the shaded light red and dark red regions highlight areas of the CMDs affected by mild and severe saturation, respectively.

4. The JWST colour-magnitude diagram of ω Cen

In the CMDs presented in Figs. 2, 3 and 4, the MS appears to be composed of two distinct components, which are clearly visible in the magnitude range from m_F150W2 ∼ 18 (m_F322W2 ∼ 17.5) down to m_F150W2 ∼ 20 (m_F322W2 ∼ 19.5). At this point, the two sequences intersect and exchange positions. The sequence that is bluer above the intersection (and redder below) corresponds to the bMS population, while the sequence that is redder above the intersection (and bluer below) corresponds to the rMS population.

Note that we have chosen to maintain the names – bMS and rMS – as they were originally introduced based on the identification of these sequences using optical filters, where the bMS appears on the blue side and the rMS on the red side of the CMD (see Bedin et al. 2004; Bellini et al. 2009). Although in the low-MS, which is the focus of our study, the two sequences appear inverted in the JWST CMD (with the rMS on the blue side and the bMS on the red side), we retain this nomenclature for consistency with previous works.

To illustrate this inversion, in Fig. 5, we have isolated a sample of bMS and rMS stars from the CMD based on optical HST photometry and shown their position on the JWST CMD. As clearly shown, the bMS crosses the rMS around m_F150W2 ∼ 20, becoming redder at fainter magnitudes. Below this crossing point, rMS stars span a broad colour range, but their distribution appears to favour bluer colours.

Fig. 5. Crossing of the bMS and rMS stars in the JWST photometry-based CMD. (a) In the HST mF814W versus mF606W − mF814W CMD, we isolated a sample of bMS (black points within blue circles) and rMS (black points within red circles) stars; the two sequences are clearly distinguishable. (b) The same samples highlighted in the JWST mF150W2 versus mF150W2 − mF322W2 CMD, using the same colour code. The shaded light red and dark red regions highlight areas of the CMDs affected by mild and severe saturation, respectively.

We note that, despite the onset of saturation at m_F150W2 ∼ 20.7, the two MSs remain clearly distinguishable in the JWST CMD, with a colour separation that appears well defined throughout the saturated regime. Although saturation introduces photometric discontinuities, it affects both sequences in a similar way and does not prevent a qualitative assessment of their morphology in this region.

Beyond this intersection (m_F150W2 > 20 or m_F322W2 > 19.5), the two main groups of stars remain distinguishable down to their respective terminations (see Figs. 2, 3 and 4). The rMS seems to terminate at m_F150W2 ∼ 24.75 (m_F322W2 ∼ 24), while the bMS ends ∼0.5 magnitude brighter (m_F150W2 ∼ 24.25 or m_F322W2 ∼ 23.5). To determine whether these terminations are real and not influenced by completeness effects, we conducted AS tests to assess the completeness of our catalogue, which will be discussed in the next section.

The scatter in photometric colour is driven by the presence of mPOPs, and it encodes the chemical spread in ω Cen. In order to interpret this information, we compared the observed CMD to theoretical isochrones. At this stage, we restricted our analysis to the overall range of m_F150W2 − m_F322W2 colour in the CMD as a function of magnitude, and ignored all finer details of the underlying colour distribution. A likelihood-based analysis of star densities in the colour-magnitude space is deferred to a dedicated future study.

At each magnitude, we defined the observed colour range as the interval between the bluest and reddest colour, for which the star density in the colour-magnitude space exceeds a certain threshold that was adjusted by hand to capture the bulk of the colour distribution and exclude outliers. The calculated range is highlighted with yellow shading in panel (a) of Fig. 6. As the basis of our comparison, we adopted two isochrones from Gerasimov et al. (2022) that the authors refer to as ‘nominal’ (NOM) and ‘High Metal High Alpha’ (HMHA). Both isochrones were calculated for the helium-enriched sub-population of the cluster with Y = 0.4 (King et al. 2012) that was first recognised in Bedin et al. (2004). The NOM isochrone uses the modal light element (C, N and O) abundances and metallicity ([Fe/H]= − 1.7) from the spectroscopic observations of giant members in Marino et al. (2012). The HMHA isochrone was derived from NOM, but with metallicity and α-enhancement (including O) tuned to match the ridgelines of the NIR and optical CMDs from HST photometry. In particular, the HMHA model adopts [M/H]= − 1.4, [α/M]= + 0.6, as detailed in Gerasimov et al. (2022). The overall range of [O/H] covered by these two isochrones spans 1 dex, and overlaps with ∼80% of the spectroscopically inferred distribution⁴.

Fig. 6. Comparison between the JWST NIR CMD of ω Cen presented in this study, and model isochrones. For clarity, we do not plot isochrones on the observed CMD directly. Instead, we compare the isochrones to the region of the colour-magnitude space occupied by the cluster, which is shaded in yellow in all panels. Panel (a) demonstrates how this region was chosen in relation to the observed CMD. Panel (b) over-plots two of the isochrones from Gerasimov et al. (2022) with both their original helium mass fraction (Y = 0.4, dashed lines) and with modified evolutionary models that incorporate a near-solar helium abundance (Y = 0.25, solid lines). Panels (c) and (d) show the Y = 0.25 NOM isochrone as well as its lower-MS variations with various chemical offsets. In the legend, ΔC, ΔO and ΔM refers to the offsets in [C/H], [O/H] and [M/H] from the chemical composition adopted for NOM in Gerasimov et al. (2022), where M only refers to metals heavier than oxygen. The shaded light red and dark red regions highlight areas of the CMDs affected by mild and severe saturation, respectively.

In order to explore the full range of helium mass fraction in ω Cen, we calculated alternative sets of evolutionary models for both NOM and HMHA with Y = 0.25 using the MESA code (Paxton et al. 2011, 2013, 2015, 2018, 2019), and following the method from Gerasimov et al. (2024). We note that the modified versions of the isochrones retain the original atmosphere models with Y = 0.4. Therefore, they have inconsistent helium mass fractions between stellar atmospheres and interiors. We anticipate the error due to this inconsistency to be insignificant since the shape of the lower MS is insensitive to Y (Libralato et al. 2024), and the effect of Y on the upper MS largely originates in the stellar interiors rather than atmospheres (Salaris & Cassisi 2005). We transformed the isochrones to the observed plane using the same parameters and procedure as in Gerasimov et al. (2022), but with a lower interstellar reddening (E(B − V) = 0.12, which is more consistent with other measurements in the literature, Harris 1996; Calamida et al. 2005) and a more up-to-date reddening law from Gordon et al. (2023).

All four model isochrones (the original and modified versions of NOM and HMHA) are compared to the observed colour spread in panel (b) of Fig. 6. The isochrones are also provided in Tables B.1, B.2, B.3, and B.4 in Appendix B. Both the isochrones and the observed CMD display a clear change of slope near m_F322W2 ≈ 19.5. This feature is known as the MS knee (Bono et al. 2010; Saracino et al. 2018), and it is caused by the reduction in the adiabatic gradient due to dissociation of molecular hydrogen (Saumon et al. 1995; Cassisi 2011). Above the knee, the effect of atmospheric chemistry on the shape of the isochrone is subtle, since high effective temperatures (T_eff ≳ 4500 K) suppress the abundances of key infrared absorbers such as H₂O. The small residual effect is driven by two primary factors. First, the mean opacity of the atmosphere reduces T_eff at higher metallicities and makes the star redder. This effect is captured in the atmosphere-interior coupling scheme adopted in Gerasimov et al. (2022). The second factor is related to the hydrogen anion (H⁻) absorption that has a minimum near 1.6 Mm, producing the characteristic ’H⁻ bump’ in the spectrum (John 1988; Sawicki 2002). The bump becomes more pronounced at higher metallicities, making the m_F150W2 − m_F322W2 colour of the star bluer.

The m_F150W2 − m_F322W2 CMD therefore provides a unique opportunity to investigate stellar interiors, since these two already small colour-metallicity dependencies due to stellar atmospheres suppress one another. For this reason, the colour of stars above the MS knee is predominantly determined by the stellar interior, and not the atmosphere. In particular, the observed colour of the star is expected to be largely determined by the mean molecular weight of the interior, which is most sensitive to the helium mass fraction, Y. Panel (b) of Fig. 6 clearly shows that Y has the dominant impact on the m_F150W2 − m_F322W2 colour, despite the NOM and HMHA isochrones having vastly different chemistries in the atmospheres. The scatter in Y of ∼0.15 is required to fully capture the observed width of the CMD above the MS knee.

The scatter in the CMD colour is noticeably wider below the MS knee due to the increasing effect of atmospheric chemistry including, most notably, infrared absorption bands of H₂O that are regulated directly by [O/Fe], and indirectly by [C/Fe] through the carbon-oxygen chemical reaction network of the atmosphere (Madhusudhan 2012). Both [C/Fe] and [O/Fe] are expected to display large member-to-member variations in GCs due to mPOPs. To determine if these variations alone are sufficient to explain the observed CMD below the MS knee, we adopted the T_eff − log(g) and T_eff-luminosity relationships from the Y = 0.25 NOM isochrone, and calculated additional model atmospheres for various values of [C/Fe]/[O/Fe] at 3500 K ≤ T_eff ≤ 5000 K in 100 K steps. The new models were calculated using the BasicATLAS/ATLAS-9/SYNTHE setup (Larkin et al. 2023; Kurucz 1970, 2005, 2014; Kurucz & Avrett 1981).

The upper temperature limit for the new models was chosen to approximately match the severe saturation limit of our photometry (also shown in Fig. 6 with dark red shading). The lower temperature limit was chosen out of practical considerations, as ATLAS atmospheres are expected to become increasingly unreliable at T_eff ≲ 3500 K (Plez 2011). Due to the low metallicity of ω Cen, the effective temperature of 3500 K corresponds to comparatively low stellar masses, varying between 0.10 M_⊙ and 0.15 M_⊙ among the four isochrones introduced above. We therefore expect the temperature range of new model atmospheres to be sufficiently wide to accommodate nearly all of our data.

For each set of chemical offsets from the composition of the NOM isochrone, we only calculated the updated synthetic spectra, and ignored the effect of altered chemistry on the atmospheric structures, atmosphere-interior coupling and stellar evolution. While these effects are expected to be subdominant on the lower MS, they may not be insignificant. For this reason, this analysis must be considered preliminary and taken with caution.

The synthetic colours and magnitudes of the new models are plotted alongside the NOM Y = 0.25 isochrone in panel (c) of Fig. 6. The figure shows that a spread in [C/Fe]/[O/Fe] alone is sufficient to explain the full width of the observed CMD; however, the red tail of the colour distribution requires a significant fraction of the stars to have simultaneously very low oxygen and very high carbon abundances, which is unlikely due to the strong carbon-oxygen correlation that is expected in GCs in general (see review in Bastian & Lardo 2018), and has been spectroscopically confirmed in ω Cen in particular (Marino et al. 2012).

Alternatively, the red tail of the colour distribution may be due to member-to-member variations in metallicity that have been spectroscopically confirmed in this GC (Marino et al. 2011). Panel (d) of Fig. 6 demonstrates that the red tail of the CMD can be approximated either by models with large but plausible metallicity offsets from NOM, or by more modest metallicity offsets in combination with adjusted light element abundances. The effect of individual metals on the m_F150W2 − m_F322W2 colour of the star at T_eff = 4000 K is shown in Fig. 7. With the exception of oxygen that directly influences infrared absorption features in the spectrum, the most important metal abundances are those of key electron donors (e.g. Mg, Ca, Na) that alter the slope of the stellar continuum.

Fig. 7. Effect of offsets in the abundances of selected metals from the chemical composition of the NOM isochrone (panels (b)–(d) of Fig. 6) on the observed mF150W2 − mF322W2 colour of the star at Teff = 4000 K. The effects are shown at three different metallicity offsets, where M only includes metals heavier than oxygen. Positive values indicate that the star gets redder when the abundance of the element is increased. The 0 (no-effect) horizontal line is shown for reference.

As seen in both panels (c) and (d) of Fig. 6, appropriately chosen isochrones can be used to approximate either tail of the distribution across most of the observed lower MS. We therefore conclude that our observations of ω Cen do not suggest that the chemical spread in the cluster varies with stellar mass.

5. Artificial stars

We used ASs to estimate the completeness of our catalogue for each of the two sequences. We first evaluated the completeness of the JWST data by generating ASs uniformly distributed across the JWST FOV. This allowed us to evaluate the completeness level for the CMD that includes all stars within the JWST field (see Fig. 2). The magnitudes in F322W2 were uniformly distributed within the range 18 < m_F322W2 < 28.5, while the corresponding F150W2 magnitudes were assigned based on fiducial lines manually defined in the m_F322W2 versus m_F150W2 − m_F322W2 CMD. One fiducial line was drawn for each of the two sequences, as shown in Fig. 8. We generated 50 000 ASs for each sequence, for a total of 100 000 ASs. The ASs were generated, detected, and measured using KS2, following the same procedures used for the real stars. An AS was considered recovered if the difference between the input and output positions was less than 1 pixel, the difference between the input and output magnitudes was within 0.75 mag (equivalent to ∼2.5log2), and it passed the same selection criteria applied to the real stars.

Fig. 8. mF322W2 versus mF150W2 − mF322W2 CMDs of our selected sample of stars, focusing on the upper-MS (panel a) and the entire MS (panel b), with the two fiducials used to generate the ASs. The dashed grey rectangle in panel b highlights the CMD region displayed in panel a. The blue fiducial corresponds to the bMS, while the red fiducial corresponds to the rMS. In the upper part of the MS, the bMS occupies the blue side of the CMD, while the rMS is on the red side. Below the intersection point, the two sequences switch positions. The shaded light red and dark red regions highlight areas of the CMDs affected by mild and severe saturation, respectively.

Panel (a) of Fig. 9 shows the injected ASs for the two sequences in the m_F322W2 versus m_F150W2 − m_F322W2 CMD, with injected stars in red and blue, and recovered stars in black. For comparison, panel (b) presents the same CMD for the real stars, with the two fiducial lines overplotted, using the same scale as in panel (a). As clearly visible from panel (a), the recovered ASs exhibit a slight asymmetry in their colour distribution, with redward tails extending from both sequences. This behaviour stems from a known effect of crowding in dense stellar fields. In particular, ASs that fall near bright real stars can be recovered as slightly brighter and redder. This systematic effect produces a characteristic skewness in the colour distribution, manifesting as redward tails.

Fig. 9. (a) Injected ASs for the two sequences (blue and red points) in the mF150W2 versus mF150W2 − mF322W2 CMD. Recovered ASs are represented in black. (b) Same as panel (a), but the black points now represent the real stars. (c)–(d) Completeness as a function of the mF150W2 and mF322W2 magnitude, respectively. The completeness for each sequence is represented using the same colour scheme as in panel (a). Shaded light red and dark red areas qualitatively mark magnitude ranges where MS stars are mildly or severely saturated in at least one filter.

The completeness, defined as the ratio of recovered stars to the total number of injected stars, is shown as a function of m_F150W2 and m_F322W2 magnitudes in panels (b) and (c), for both sequences. In these panels, the completeness for each sequence is represented using the same colour scheme as in panel (a). In this and the other figures of this section, shaded light red and dark red areas qualitatively mark magnitude ranges where MS stars are mildly or severely saturated in at least one filter. Completeness estimates in these magnitude ranges should be interpreted with caution.

From panel (a), it is evident that the recovered ASs extend well beyond the termination points observed in the real stars. The recovered ASs reach the same depth for both sequences, showing no signs of the differing termination magnitudes observed in the real stars. The completeness levels for the two sequences, shown in panels (b) and (c), are nearly identical. At the observed termination magnitudes of the sequences in the real stars (m_F322W2 ∼ 24 for the rMS and m_F322W2 ∼ 23.5 for the bMS), the completeness remains relatively high, around ∼40%.

These results confirm that the observed terminations in the two sequences are real and not influenced by completeness limitations. The absence of a significant difference in completeness levels between the two sequences further supports the conclusion that the earlier termination of the bMS is an intrinsic feature rather than a consequence of completeness effects.

To evaluate the completeness level for the CMD that includes all stars in the overlapping FOV of both HST and JWST (see Fig. 3), we generated an additional set of 100 000 ASs (50 000 for each sequence), uniformly distributed across the common area. The resulting completeness levels are shown in Fig. 10 as a function of the m_F150W2 (panel a) and m_F322W2 (panel b) magnitudes. The completeness values are nearly identical to those obtained when considering the entire JWST field.

Fig. 10. Completeness as a function of mF150W2 (panel a) and mF322W2 (panel b) magnitudes for the two sequences in the overlapping region of the HST and JWST fields.

To incorporate PM selection into our completeness evaluation and assess the completeness for the stars shown in panel (e) of Fig. 4, we performed the ASs procedure for the HST data as well. We used the same positions and F150W2 and F322W2 magnitudes defined above and determined the corresponding F606W and F814W magnitudes using fiducial lines manually defined in the m_F322W2 versus m_F814W − m_F322W2 and m_F322W2 versus m_F606W − m_F322W CMDs. We ran the ASs stars separately for the JWST and HST datasets and then assessed the displacement of the stars between the two epochs.

Since ASs are generated with identical positions in both epochs, any observed displacement is solely due to spurious offsets caused by noise (e.g. uncertainties in PSF modelling, coordinate transformations, cosmic ray hits, detector cosmetics, etc.), which introduces slight shifts in the recovered positions across epochs. To be considered recovered, an AS must also satisfy the same PM selection criteria applied to real stars, as shown in panel (d) of Fig. 4. By applying this PM selection to ASs, we estimated the fraction of real stars excluded due to these positional offsets, allowing us to incorporate this effect into the overall completeness calculation.

The completeness obtained after applying the PM selection is shown in Fig. 11. As expected, the completeness of both sequences is lower and declines more rapidly towards fainter magnitudes compared to Fig. 10, reaching zero at a brighter magnitude. However, at the termination points of the sequences for real stars (m_F322W2 ∼ 24 for the rMS and m_F322W2 ∼ 23.5 for the bMS; see panel (d) of Fig. 4), the completeness remains relatively high, around ∼30–40%, confirming the reliability of the sequences termination. Given the consistent completeness between the two sequences, we adopt the completeness values obtained for the rMS in the following analysis.

Fig. 11. Same as Fig. 10 but including the PMs selection.

6. Luminosity and mass function

In this section, we present the LF and MF of ω Cen and its sequences using JWST photometry.

6.1. Combined luminosity and mass functions

We derived the combined (i.e. without matching individual stars to sub-populations) LF and MF of ω Cen by fitting a forward model to the observed m_F150W2 and m_F322W2 magnitudes of confirmed members. The likelihood of compatibility between the forward model and a set of photometric measurements with uncertainties, {m_i, σ_i}, is given by the following equation:

$$\begin{array}{c}\hfill \mathcal{L}={\displaystyle \prod _i}{\displaystyle {\int }_{m_{\min }}^{\infty }}\varphi (m_i)\mathrm{comp}(m_i)Z(x;m_i,{\sigma }_i)dx,\end{array}$$(1)

where ϕ(m) is the LF (i.e. the likelihood of observing a member with magnitude m), comp(m) is the photometric completeness at magnitude m, and Z(x; m, σ) is the probability of measuring magnitude x, given the true magnitude m and the photometric error σ. We took Z(x; m, σ) to be a normal distribution with the mean of m and the standard deviation of σ. The brightest considered magnitude, m_min, was set to the saturation limit of the instrument, which we conservatively set to m_min = 19 for both F150W2 and F322W2.

We considered five different families of LF models (ϕ(m)) listed below:

1. Single population with single-component power-law MF. To construct this model, we used the mass-luminosity relationship (MLR) from only one of the theoretical isochrones shown in panel (b) of Fig. 6. We took the MF to obey a single-component power-law distribution as follows:

   $$\begin{array}{c}\hfill \xi (M)\propto M^{-{\alpha }_{\mathrm{h}}},\end{array}$$(2)

   where ξ(M) is the MF (i.e. the likelihood of observing a member with mass M), and α_h is the slope of the power law. ξ(M) was converted to ϕ(m) using the MLR of the adopted isochrone as follows:

   $$\begin{array}{c}\hfill \varphi (m)\propto \xi (M(m))\left|\frac{dM(m)}{\mathit{dm}}\right|,\end{array}$$(3)

   where M(m) is the MLR, and the derivative was evaluated using a linear spline interpolator.

2. Single population with two-component broken power-law MF. As before, only one isochrone is used to derive the LF, but the MF is allowed to have a break in the power law:

   $$\begin{array}{c}\hfill \xi (M)\propto \{\begin{array}{cc}{M}^{-{\alpha }_{\mathrm{h}}},\hfill & \text{if}M>M_{\mathrm{pb}}\hfill \\ M^{-{\alpha }_{\mathrm{l}}},\hfill & \text{if}M\le M_{\mathrm{pb}}\hfill \end{array},\end{array}$$(4)

   where α_h and α_l are the power law slopes above and below the break, respectively, and M_pb is the stellar mass at which the break occurs.

3. Mixture of two populations with identical single-component power-law MFs. This family of models utilises the single-component MF in Eq. (2); however, the LF is composed of two populations with MLRs adopted from two distinct isochrones in panel (b) of Fig. 6. The contributions of the two populations to ϕ(m) are added with weights given by γ and (1 − γ), where 0 ≤ γ ≤ 1 is the mixing ratio. Both populations are assumed to have identical MFs.

4. Mixture of two populations with identical two-component broken power-law MFs. As before, but the shared MF of the two populations has a power law break (Eq. 4).

5. Mixture of two populations with distinct single-component power-law MFs. These models consider two populations with single-component MFs in Eq. 2; however, the slope of the power law is allowed to vary between the populations. We denote the slope of the MF of population 1 as α_h, and the slope of population 2 as α₂.

We did not consider more complicated LF models (e.g. two-population models with distinct broken power-law MFs or three-population models) to avoid over-fitting the observations; however, three-population configurations and arbitrary functional forms of the MF are explored in Section 6.2, where the observed CMD sequences are analysed separately.

The logarithmic likelihood in Eq. 1 was maximised using the Goodman-Weare (Goodman & Weare 2010) Markov chain Monte Carlo (MCMC) sampling with respect to the free parameters defined above (α_h, α_l, α₂, M_pb and γ, as appropriate for each family of LF models). We used 32 walkers and 3000 accepted steps per walker. The first 10% of the accepted steps were discarded as burn-in. For single-population LF models, we ran separate MCMC chains for each of the four isochrones in panel (b) of Fig. 6. For LF models with two populations, we considered every possible pair of isochrones, except the ones where both populations have Y = 0.4, as such configurations are clearly unphysical. Furthermore, we ran separate chains for the LFs in F150W2 and F322W2 bands to verify the internal consistency of our results. The best-fit values of the free parameters were taken as the medians of the corresponding MCMC posteriors. The upper and lower 1-sigma asymmetric errors were taken as half-differences between the median and the 97.7 and 2.3 percentiles. The results are summarised in Table 1.

We note that our model fitting method does not require binning of data, and is free of issues associated with histogram-based approaches (e.g. Sollima et al. 2007b; Bedin et al. 2024b; Gerasimov et al. 2024) such as the dependence of results on the bin size and incorrect bin placements due to photometric errors. However, it is possible to visualise the best-fit model in the histogram form by integrating the completeness-corrected LF (ϕ(m_i) comp(m_i) in Eq. (1)) in the chosen magnitude bins. This has been carried out in 15 uniform magnitude bins to produce panels (a) and (b) of Figure 12. These panels show the observed member counts in each bin, as well as the best-fitting (based on the maximum likelihood in the MCMC chain for F150W2) LF models from each of the five families described above. In all cases, the best-fitting models are the ones that use the Y = 0.25 NOM isochrone (for single-population models) or both Y = 0.25 and Y = 0.4 NOM isochrones (for mixed population models). The errors shown in Figure 12 were estimated from the scatter among a large number of synthetic datasets that were generated using the best-fit LF model and the estimated completeness and photometric errors of our observations.

We also derived the goodness-of-fit (GoF) for each model as the average absolute difference between the observed and modelled histogram counts, normalised by the estimated count errors. The GoF of each considered LF model is shown in Table 1. Note that for 15 magnitude bins, assuming Gaussian random errors and no systematic errors, the expected GoF is ≈0.80 ± 0.16⁵. Models with GoFs significantly higher than this value are likely over-fitted, whereas those with lower values are likely under-fitted.

The key conclusions that can be drawn from Table 1 and Fig. 12 are summarised below:

Fig. 12. Observed LFs and derived MF of ω Cen. Panels (a) and (b) show the observed LFs of ω Cen in F150W2 and F322W2 bands, respectively. The LFs are represented as histograms with 15 uniform bins. Over-plotted in both panels are the best-fitting LF models (taken to be the model with the smallest GoF in F150W2) from each of the five families listed in Table 1. The model curves include completeness effects. Panel (c) shows the MF of the cluster based on the inferred stellar masses from the observed F150W2 magnitudes and the Y = 0.25 NOM isochrone, shown in panel (b) of Fig. 6. Unlike panels (a) and (b), the inferred counts in panel (c) have been corrected for photometric completeness. Over-plotted are the best-fitting single-component and two-component power-law MFs. In the legends, pop stands for population, comp for component, and PL for power law. The light red areas qualitatively mark magnitude ranges where MS stars are mildly saturated in at least one filter. The LF and MF values in these regions should be interpreted with caution.

1. Out of single-population models with unbroken single-component MFs, NOM25 and NOM40 provide the best fits to the data (smallest GoF) in F150W2 and F322W2, respectively. The latter case is clearly unphysical, since the majority of members in the observed field are expected to belong to rMS that has a near-solar Y ≈ 0.25. As we argue later in this section, this contradiction most likely arises because ω Cen does have a break in the MF within the observed mass range (∼0.1 − 0.5 M_⊙), and because super-solar values of Y mimic such a break in the LF (higher helium mass fractions create a discontinuity in the MLR near 0.3 M_⊙, where stars become fully convective). If we only consider single-population unbroken-MF models with Y = 0.25, then NOM25 provides the best fit in both F150W2 and F322W2. The inferred MF slopes in both bands (∼0.7 − 0.8) are consistent within uncertainties and are also consistent with similar analyses in the literature (e.g. Sollima et al. 2007b; Gerasimov et al. 2022).

2. LF models with broken power-law MFs provide better fits to the data than the unbroken-MF models for all considered configurations. In most cases, adding a break to the MF improves the GoF by a factor of ∼2, and brings it in line with the expected value of ≈0.80 ± 0.16 for a well-fitting model. Out of single-population broken-MF models, NOM25 offers the best fit in both F150W2 and F322W2. In both bands, a power law break around M_pb ≈ 0.2 M_⊙ is observed with a steep bottom-heavy MF above the break (α_h > 1 at M > M_pb) and a nearly flat MF below the break (α_l < 0.2 at M ≤ M_pb). Panel (c) of Fig. 12 shows a histogram of the completeness-corrected inferred stellar masses using the NOM25 broken-MF model. The break in the power law is clearly seen.

3. As is emphasised above, the apparent break in the MF may actually be the result of a discontinuity in the MLR due to the transition from partially to fully convective stellar interiors. In our models, this transition occurs near 0.32 M_⊙ for Y = 0.25, and 0.27 M_⊙ for Y = 0.4. We however note that the best-fit M_pb for the broken-MF NOM25 model was estimated as 0.19 ± 0.03 M_⊙, i.e. over three sigma below the transition. This further corroborates that the observed MF break is likely genuine. On the other hand, all broken-MF single-population models with Y = 0.4 place the break near 0.3 M_⊙. For these models, the inferred break is an artefact caused by the discontinuity in the MLR due to the transition from partially to fully convective stellar interiors.

4. While introducing a break in the MF improves the GoF considerably, adding a second population to the mixture has a very small effect. This is apparent in panels (a) and (b) of Fig. 12, where the broken LF models with single and mixed populations are nearly indistinguishable from each other.

5. The best-fitting mixed population models are the combinations of NOM25+NOM40 and NOM25+HMHA25 populations with a common two-component broken MF. In the former case, the mixing ratio was estimated as 0.7 ± 0.1, suggesting that approximately 20 − 40% of cluster members are expected to be significantly helium-enriched (bMS). This fraction is consistent with the fraction of observed bMS at this distance from the centre of ω Cen (Scalco et al. 2024a; Bellini et al. 2009). This result is a significant improvement over a similar analysis of HST observations in Gerasimov et al. (2022), where only an upper limit on the population mixing ratio could be derived. In the following sub-section, we verify the calculated mixing ratio by analysing the individual sequences of ω Cen separately.

6. The best-fitting broken-MF models in both single-population and mixed configurations reach the target GoF for a well-fitting model (≈0.80 ± 0.16). Addition of further degrees of freedom to the model would result in over-fitting. For this reason, we do not consider mixed models with distinct two-component MFs, or mixed models with more than two populations. A far larger sample size is required to explore these more sophisticated models.

7. Allowing a two-population LF model to have distinct MFs in each population has a very small effect on the GoF, compared to the corresponding LF models with identical single-component MFs. For nearly all distinct-MF models, the mixing ratio was found to be consistent with γ = 1 within 2-sigma bounds, and no reliable constraints on α₂ could be derived. A notable exception is the NOM25+HMHA40 mixture, for which all parameters are well constrained, and the best-fit values are consistent between F150W2 and F322W2. In this model, both populations were found to have similar MF slopes (α_h ≈ α₂ ≈ 0.5), which are comparable to the MF slopes calculated for single-population models with Y = 0.4. It is therefore likely that the best-fit parameters of the NOM25+HMHA40 model are an artefact of the discontinuity in the MLR at the transition between fully and partially convective stellar interiors.

Fig. 13. Procedure for estimating the LFs of the sequences in ω Cen in the low-MS. Panel (a) shows the mF322W2 versus mF150W2 − mF322W2 CMD, focusing on the region where the separation between the components is most evident. The red and blue lines indicate the fiducials of the rMS and bMS, respectively, and are used to construct the verticalised CMD shown in panel (b). Panels (1,2c) through (1,2l) display the histogram distribution of ΔmF150W2 − mF322W2 for stars across eight magnitude intervals, as defined by the grey lines in panels (a) and (b). In panels (1c) to (1l), the red and blue regions corresponding to the two sequences are highlighted. These regions are used to estimate the LFs presented in Fig. 14, with values reported in Table 2. Panels (2c) to (2l) present the same histograms with the best-fitting multi-Gaussian models overlaid: three-Gaussian fits in panels (2c) to (2i) and a two-Gaussian fit in panel (2l). The individual Gaussian components are shown in red, green, and blue, while the combined fit is plotted in grey. The obtained LFs are shown in Fig. 15 and the values are reported in Table 3.

Fig. 14. (a) LF derived using the regions defined in panels (1c) to (1l) of Fig. 13. (b) Population ratio of stars associated with the two sequences. In both panels, the red-shaded region indicates the saturation limit, while the grey-shaded region marks the threshold beyond which the three populations can no longer be distinguished.

Fig. 15. (a) LF obtained from the three-Gaussian fit shown in panels (2c) to (2l) of Fig. 13, highlighting the drop of the bMS in the final bin (mF150W2 ∼ 23.75). (b) Population ratio of stars associated with the three sequences. As in Fig. 14, both panels feature a red-shaded region indicating the saturation limit and a grey-shaded region marking the threshold beyond which the three populations can no longer be distinguished.

6.2. Luminosity and mass functions for individual stellar populations

In Section 6.1, we analysed the combined LF and MF of ω Cen (where ‘combined’ refers to the local LF and MF obtained by merging all subpopulations). However, this analysis did not fully exploit all the information available in the CMD, particularly the colour of each star relative to the overall colour distribution, which provides an indication of its likely population membership. As a result, we were limited in our ability to explore more complex models (such as two-population scenarios with distinct broken power-law MFs, or three-population models) and to examine in detail the differences in MF slopes among the individual populations. To overcome these limitations, we adopted a more direct approach, focusing on the LF and MF of the individual stellar populations over a narrower mass range, where the sequences are partially distinguishable. Unlike the method used in Section 6.1, here we adopted a binned approach, as it allows for a more accurate separation of the individual sequences and facilitates a clearer comparison of their respective LFs.

We derived the LFs along the low-MS of the two populations, applying a methodology similar to that outlined in Scalco et al. (2024d). The method is illustrated in Fig. 13. Panel (a) of Fig. 13 shows the m_F322W2 versus m_F150W2 − m_F322W2 CMD for our selected sample of stars, in the magnitude range 20 < m_F322W2 < 24, where the separation of the two components is most evident. Note that although completeness remains reliable down to m_F322W2 ∼ 26, the statistical limitation of the sample prevents us from classifying or distinguishing additional sequences beyond m_F322W2 ∼ 24 (see Figs. 2, 3, 4). We utilised the previously defined fiducials (see Fig. 5) for the two sequences (shown in red and blue in panel (a) of Fig. 13) to verticalise the diagram. For each star, we computed the verticalised colour as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{F}150\mathrm{W}2-\mathrm{F}322\mathrm{W}2}=\frac{X-X_{\mathrm{red}\mathrm{fiducial}}}{X_{\mathrm{blue}\mathrm{fiducial}}-X_{\mathrm{red}\mathrm{fiducial}}},\end{array}$$(5)

where X = m_F150W2 − m_F322W2 and X_{red fiducial}, X_{blue fiducial} are are the colours of the red and blue fiducials, respectively, evaluated at the star’s m_F322W2 magnitude. The resulting verticalised diagram is shown in panel (b) of Fig. 13. We defined 8 magnitude bins in m_F322W2, each 0.5 magnitudes wide (indicated by the grey horizontal lines in panels (a) and (b)). The histogram of the verticalised colour for each magnitude bin is shown in panels (1,2c) to (1,2l). We defined two regions: one with −0.4 < Δ_{mF150W2 − mF322W2} < 0.6 corresponding to the rMS, and another with 0.6 < Δ_{mF150W2 − mF322W2} < 1.4 corresponding to the bMS. These regions are represented in panels (1c) to (1l) in red and blue, respectively. For each bin, we counted the number of stars within each of the two defined regions. The resulting values, corrected for completeness, are shown in panel (a) of Fig. 14, using the same colour scheme as in Fig. 13, with error bars representing Poisson errors (this convention for the error bars is adopted throughout this figure and in all subsequent LF and MF plots). As observed, the two LFs exhibit a similar shape, with the rMS containing a larger number of stars in each magnitude interval. The bMS LF declines more rapidly than the rMS LF at fainter magnitudes (m_F322W2 > 22.5) due to the earlier termination of the bMS sequence.

Panel (b) of Fig. 14 presents the population ratio between the two sequences. The error bars represent uncertainties derived through standard error propagation – this approach is also adopted in all subsequent figures displaying population ratios. We calculated the weighted mean of the ratios, finding that the rMS accounts for 75%±1%, while the bMS accounts for 25%±2% of the total number of stars. These values are consistent with those previously reported in literature (Scalco et al. 2024a; Bellini et al. 2009, bMS/rMS ∼ 0.33, for a radial distance corresponding to the field analysed in this study, ∼2 − 3r_h), and agree with the results discussed in Section 6.1. The data presented in Fig. 14 are listed in Table 2.

To estimate the LFs of the two sequences while accounting for possible contamination between them, we used an alternative method, detailed in Section 4.2 of Scalco et al. (2024d). This method involves determining the number of stars in each bin and in each sequence by fitting the histograms with a multi-Gaussian model, where the number of Gaussians corresponds to the number of sequences. In our case, we initially aimed to fit the histograms with two Gaussian components, one for each sequence. However, upon closer inspection of the verticalised diagram in panel (b) and the histograms in panels (1c, 2c) to (1l, 2l), we observed evidence of a possible third component situated between the two MSs. As discussed in the introduction, ω Cen hosts a highly complex system of mPOPs. While the bMS and rMS constitute the dominant populations, additional stellar groups have been identified, each exhibiting substructures (including the rMS and bMS), resulting in at least 15 distinct populations within ω Cen (Bellini et al. 2017a; Scalco et al. 2024a). The MS presented in panels (a) and (b), and the relative histograms shown in panels (1,2c) to (1,2l), are the result of the contribution of all these stellar populations, which overlap and blend with each other. This complexity necessitates a more precise method to determine the optimal number of Gaussian components to fit.

To address this, we applied a Gaussian mixture model (GMM) to the Δ_{mF150W2 − mF322W2} distribution of stars in each bin, using the expectation-maximisation algorithm from the scikit-learn package (Pedregosa et al. 2011). To find the best model, we calculated the Bayesian Information Criterion (BIC) for models with one to five Gaussian components. A three-Gaussian model provided the lowest BIC for most bins, except for the last two bins (those with, respectively, 23 < m_F322W2 < 23.5 and 23.5 < m_F322W2 < 24) where a two-Gaussian and a one-Gaussian model were optimal, respectively. In these final bins, the rightmost sequence is barely detectable in the penultimate bin and entirely absent in the last, leaving only the leftmost sequence in the last bin, along with a faint remnant of the middle sequence. Based on these results, we adopted a three-Gaussian model for all bins except the last one, where a two-Gaussian model was used. The GMM fits are shown in grey in panels (2c) to (2l), with the individual components represented in red, green, and blue. Hereafter, we refer to the leftmost and rightmost components as the rMS and bMS, respectively, while the middle sequence are referred to as the green MS (gMS).

For each bin, the number of stars in each sequence was estimated from the area under the corresponding Gaussian component. The resulting LFs, corrected for completeness, are illustrated in panel (a) of Fig. 15. As shown, the three LFs follow a similar shape. The rMS is the most populous sequence, the bMS is the least populous, and the gMS lies in between. The only exception occurs at m_F322W2 = 20.75, where a slight peak in the bMS LF temporarily exceeds that of gMS. In the lower part, the LF for the bMS shows a drop, reaching zero in the final bin as the sequence seems to end. Panel (b) shows the population ratios of the three sequences across each magnitude bin. We calculated the weighted mean of the ratios, finding that the rMS accounts for 52%±2%, the gMS for 26%±1%, and the bMS for 22%±2% of the total number of stars. The data presented in Fig. 15 are listed in Table 3.

Fig. 16. (a) MFs of the three identified sequences, along with interpolating lines in the mass range 0.1 < M < 0.3 M⊙ (solid lines). The corresponding power-law slopes are indicated in the plot. For the gMS and bMS, we also show fits in the narrower mass range 0.2 < M < 0.3 M⊙ (dashed lines), with their respective slopes also reported. The black histogram represents the combined MF, as shown in panel (c) of Fig. 12. The solid black line and dashed yellow line correspond to the single-component and broken (two-component) power-law fits to the combined MF, respectively, with their slopes provided for reference. (b) Population ratio of stars associated with the three sequences as a function of stellar mass. The shaded light red and dark red regions highlight areas affected by mild and severe saturation, respectively.

Table 3 also reports the dispersions (σ) of the three Gaussian components used to model the colour distribution of the MS stars. To assess whether the observed colour spread of the three components arises from intrinsic star-to-star variations in C and O abundances – rather than being solely due to photometric errors – we performed a similar analysis on the ASs shown in panel (a) of Fig. 9. We verticalised the two ASs sequences (black points in panel (a) of Fig. 9) using the same fiducials as for the real stars, and divided the verticalised CMD into 0.5-mag bins in m_F322W2, over the same magnitude range used for the real data (20 < m_F322W2 < 24). In each bin, we applied the GMM to estimate the colour dispersion of the two components. The resulting σ values are listed in Table 3. These dispersions are systematically smaller than those measured from the real data, except in the two faintest bins, where photometric errors start to become significant – particularly for the bMS. It is worth noting that the dispersions for the two injected sequences are very similar across the full magnitude range, diverging only in the last two bins. This comparison confirms that the colour dispersion observed in the real stars cannot be explained by photometric errors alone, and is consistent with intrinsic variations in chemical composition, most notably in C and O abundances.

As was pointed out by the referee, the red tail in the artificial stars discussed in Section 5 could, in principle, contribute to the apparent presence of the gMS. However, we note that such red tails involve only a small fraction of the recovered stars, which is insufficient to account for the number of stars observed in the gMS. Moreover, both the colour dispersion and the extent of these tails are significantly smaller than those measured in the real data, particularly in the region occupied by the gMS. We therefore find no indication that the gMS could be the result of this effect.

To derive the MFs of the three stellar sequences, we adopted a statistical approach. In each magnitude bin defined in Fig. 13, we randomly selected three sub-samples of stars and assigned them to the three populations. The number of stars in each sub-sample was determined according to the area under the corresponding Gaussian component. The magnitudes of the stars were then converted into masses using the isochrones shown in Fig. 6. In particular, we adopted the NOM isochrone with Y = 0.25 for stars assigned to the rMS and gMS, and the NOM isochrone with Y = 0.40 for those assigned to the bMS (see Tables B.1 and B.2 for the corresponding tabulated MLRs). We then computed the MF for each sequence by counting the number of stars in logarithmic mass bins. This entire procedure was repeated 1000 times, and the final MFs were obtained by taking the median number of stars in each mass bin across all realisations.

The resulting MFs for each sequence are shown in Fig. 16, alongside the combined MF previously presented in panel (c) of Fig. 12. The MFs of the individual populations cover a mass range of ∼0.30–0.09 M_⊙ for the rMS and gMS, and ∼0.27–0.08 M_⊙ for the bMS. Each MF was fitted with a single power-law function (see Eq. (2)) over the mass range 0.1 < M < 0.3 M_⊙. The best-fit lines and the corresponding slope values (α) are shown in the figure.

Among the populations, the gMS exhibits the steepest slope (α_g = 0.77 ± 0.12), while the rMS and bMS show similar and flatter slopes (α_r = 0.38 ± 0.04 and α_b = 0.41 ± 0.14, respectively). The MF of the rMS – the most populated sequence – is well characterised by a single power law, indicating a reliable MLR. In contrast, the MFs of the gMS and bMS are noisier and display a noticeable flattening around ∼0.2 M_⊙, which also causes the flattening observed in the combined MF at the same mass. This behaviour could indicate a real drop in the number of stars formed below ∼0.2 M_⊙, or it might reflect uncertainties in the adopted MLR, particularly at low masses. It is also interesting to notice that for masses smaller than 0.2 M_⊙, both the gMS and the bMS are characterised a MF steeper than that of the rMS. This steep slope is responsible for the steepening of the combined MF.

It is important to note that dividing the MS into only three components represents a simplification, given the high level of complexity in the stellar populations of ω Cen. Moreover, the sequence identified as the rMS in the optical HST CMD evolves into a broader colour distribution in the NIR JWST CMD at magnitudes fainter than m_F150W2 ≃ 20 (or m_F322W2 ≃ 19; see Fig. 5), with a main peak corresponding to the red Gaussian component (Fig. 13) and an extended tail that overlaps with the other two components. These features suggest a higher level of complexity in the stellar population structure, and a more refined analysis would require a more sophisticated and detailed approach.

Finally, panel (b) displays the population ratios of the three sequences across each mass bin. We computed the weighted mean of these ratios, finding that the rMS comprises 50%±1% of the total number of stars, while the gMS and bMS account for 27%±1% and 23%±1%, respectively. These values are consistent with those derived from the LFs (see Fig. 15) and listed in Table 3. The data presented in Fig. 16 are listed in Table 4.

6.3. Note on the effect of the helium mass fraction

As shown in Section 4, the stars of the bMS population are likely helium-enriched compared to their rMS counterparts with the difference in Y between the two populations reaching ∼0.15. In general, a star with enriched helium content has a higher mean molecular weight in the interior, and is more luminous than an equally massive star with a near-solar Y. However, a helium-enriched population of stars would also be able to sustain nuclear fusion at lower stellar masses, thereby lowering the hydrogen-burning limit. These two effects largely suppress each other, making the lower MS almost completely insensitive to Y.

In other words, the observed LF is determined not by the MLR, but by the slope of the MLR with respect to stellar mass (see Eq. (3)). While Y has a large impact on the MLR, its effect on the slope is subtle compared to the effect of, for example, atmospheric chemistry. For this reason, it is unlikely that the high helium content of the bMS is responsible for its early termination in the CMD compared to the rMS. This is further illustrated in Table 5, where we provide the hydrogen-burning limits and their corresponding magnitudes for the isochrones considered in this study. Note that the large effect of Y on the former does not propagate to the latter.

7. Summary

We presented the first study of the most massive GC in our Galaxy, ω Cen, using recently acquired JWST proprietary data (GO-5110; Bedin et al. 2024a). The results of our investigation are summarised as follows.

We constructed the cluster’s NIR CMDs using JWST photometry. Figure A.1 displays these CMDs across broad colour and magnitude ranges, with the onset of mild and strong saturation clearly marked (see Section 2 for further details). In these CMDs, two sequences, – corresponding to the optical bMS and rMS populations of the cluster – can be clearly separated from m_F150W2 ∼ 18 (m_F322W2 ∼ 17.5) down to m_F150W2 ∼ 20 (m_F322W2 ∼ 19.5), where they intersect and change their relative positions (see Fig. 5). Below this intersection, the bMS shifts towards redder colours, while the rMS extends over a broad colour range, predominantly favouring bluer colours. The sequences remain distinguishable down to the end of the MS. The sequence corresponding to the bMS stars appears to end – at least – ∼0.5 magnitudes brighter than the sequence associated with the rMS population. We performed AS tests to evaluate the completeness of our data and assess whether the observed sequence terminations are intrinsic or driven by incompleteness. The results clearly indicate that the earlier termination of the bMS, which occurs at significantly brighter magnitudes compared to the rMS, is an intrinsic property of the population and not an artefact of completeness effects.

Our comparison with theoretical isochrones confirms that the colour spread in the CMD is primarily driven by variations in helium (Y) abundance above the MS knee and by atmospheric chemistry below the knee. We found that a scatter in Y of approximately 0.15 is necessary to fully capture the observed width of the CMD above the MS knee, reinforcing the idea that helium enhancement plays a dominant role in shaping the stellar populations of ω Cen.

Below the MS knee, the broader colour distribution is largely influenced by variations in oxygen and carbon abundances, which regulate the strength of H₂O absorption bands in the infrared. While the observed spread can be approximated by models with a range of [C/Fe] and [O/Fe], reproducing the full extent of the red tail requires extreme chemical compositions that are unlikely given the known abundance correlations in GCs. Instead, we find that a combination of metallicity variations and light element abundance changes provides a more plausible explanation for the observed CMD morphology.

By combining our JWST data with archival HST observations, as recently presented in Scalco et al. (2024c), we calculated PMs and conducted a membership analysis. We then derived the combined LF and MF of ω Cen, independent of individual populations, by fitting a forward model to the observed magnitudes of confirmed members using an MCMC method. We considered five different families of LF models: (1) a single population with a single-component power-law MF, (2) a single population with a two-component broken power-law MF, (3) a mixture of two populations with identical single-component power-law MFs, (4) a mixture of two populations with identical two-component broken power-law MFs, and (5) a mixture of two populations with distinct single-component power-law MFs. We explored various combinations of isochrones for the populations and evaluated their contributions to the LF.

Our results indicate that the single-population broken power-law models, particularly those using the Y = 0.25 NOM isochrone, provide the best fit. The MF exhibits a break around 0.2 M_⊙, with a steep slope above the break and a flatter slope below it. Adding a second population did not substantially improve the fit. The observed LF is predominantly influenced by a single population with an average chemical composition. The mixed population models suggest that approximately 20–40% of cluster members may be helium-enriched.

We then analysed the LF and MF of individual MS populations in ω Cen by focusing on a narrower mass range (∼0.3–0.08 M_⊙) where the sequences are more clearly separated. The LFs of both sequences exhibit similar shapes, with the rMS containing a larger fraction of stars, accounting for 75%±1% of the total population. The bMS LF declines more rapidly than the rMS LF at fainter magnitudes (m_F322W2 > 22.5) due to the earlier termination of the bMS sequence. The bMS-to-rMS number ratio of ∼0.33 in this external field agrees with values previously reported by Scalco et al. (2024a), Bellini et al. (2009) for fields at a similar radial distance from the cluster centre.

To refine our LF estimates and account for potential contamination between the sequences, we adopted an alternative approach. We analysed the verticalised CMD by fitting the star distributions in different magnitude bins with a multi-Gaussian model. By closely inspecting the verticalised CMD, we identified evidence of a possible third sequence between the bMS and rMS. To determine the optimal number of Gaussian components for each magnitude bin, we applied a GMM. The GMM results favoured three Gaussian components for most bins. The LFs were then estimated based on the area under each Gaussian component.

The derived LFs show similar shapes for all three components. In the lower part of the LFs, the bMS exhibits a drop, reaching zero in the final bin. The rMS is the most populous, contributing 52%±2% of the stars, followed by the gMS at 26%±1%, and the bMS at 22%±2%.

We then estimated the MF for each of the three sequences in the mass range between 0.3 M_⊙ and 0.08 M_⊙. The rMS is well-fitted by a single power law in this mass range. The gMS and the bMS, on the other hand, are characterised by a flattening in their MFs for masses < 0.2 M_⊙ and a slope for masses > 0.2 M_⊙ that, for both populations, is steeper than that of the rMS. The steepening of the combined MF for masses smaller larger than 0.2 M_⊙ and the flattening at masses < 0.2 M_⊙ are linked to the variation in the MF of the gMS and bMS. The flattening around ∼0.2 M_⊙ for the gMS and the bMS may reflect a real drop in the number of low-mass stars or might be due to uncertainties in the adopted MLR. The gMS and bMS MF steeper than the rMS for masses > 0.2 M_⊙ may be the result of the effects of internal two-body relaxation on the radial variation in the MF for the two populations initially more centrally concentrated than the rMS (see Livernois et al. 2024, for further discussion and details; see also e.g.) for the effects of two-body relaxation on the evolution towards energy equipartition and Ziliotto et al. (2025), for an observational study of energy equipartition in ω Cen]2018MNRAS.476.2731V. We emphasise, however, that the narrow mass range and uncertainties in the MF derivation call for caution in the interpretation of the observational results. Additional studies covering a broader mass range will be necessary to further investigate the differences between the MF of multiple populations in ω Cen.

This work represents the first study of the multiple populations in ω Cen using JWST data. As is well known, ω Cen hosts up to 15 distinct stellar populations in its core. In this study, we focused on the LFs of the sequences identifiable with JWST photometry. However, the complexity of the multiple populations in this cluster requires a more in-depth analysis to fully disentangle and characterise its stellar populations. This is particularly crucial when combining JWST and HST photometry, as demonstrated in a recent study (Scalco et al. 2024b). A comprehensive investigation utilising this combined approach will be the subject of a forthcoming paper (Gerasimov et al., in preparation).

Acknowledgments

Michele Scalco and Luigi Rolly Bedin acknowledge support by MIUR under the PRIN-2017 programme #2017Z2HSMF, by INAF under the PRIN-2019 programme #10-Bedin, and by INAF under the WFAP project, f.o.:1.05.23.05.05. JA, AB, and EV acknowledge support from JWST program GO-5110. We thank our Referee, Dr. Ricardo Salinas, for his prompt and insightful review of our manuscript. His constructive comments and suggestions have significantly improved the clarity and quality of this work. This research is based on observations made with the NASA/ESA Hubble Space Telescope and the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contracts NAS 5–26555 for HST and NAS 5-03127 for JWST. These observations are associated with programs HST-GO-14118+14662 and JWST-GO-5110. The data described here may be obtained from http://doi.org/10.17909/7prx-4905.

References

Appendix A: Overview of the JWST photometric catalogue

Fig. A.1. Colour-magnitude diagrams of ω Cen using JWST filters. Panel (a) shows the mF150W2 versus mF150W2 − mF322W2 CMD, while panel (b) displays the mF322W2 versus mF150W2 − mF322W2 CMD, including all stars within the JWST FOV. In both panels, black dots represent stars that pass both the photometric quality and PM selection criteria, while grey dots indicate the remaining stars. The shaded light red and dark red regions highlight areas affected by saturation: the light red region corresponds to saturated photometry, while the dark red region marks severe saturation, where even frame zero is saturated.

Appendix B: Mass-luminosity relations

All Tables

All Figures

Fig. 1. NIRCam JWST field from the GO-5110 programme (red) and the ACS/WFC HST field from the GO-14118+14662 programme (blue) overlaid on a DSS image of ω Cen. The overlap region between the two datasets is highlighted in magenta. Units are in arcminutes from the cluster’s centre. The dashed white and green circles represent the core radius (rc = 2′⋅37) and the half-light radius (rh = 5′⋅00), respectively, with additional circles marking 2 rh, 3 rh, and 4 rh. The values of rc and rh are from Trager et al. (1995), McLaughlin & van der Marel (2005) as reported in the Harris (1996, 2010) catalogue.

In the text

Fig. 2. Colour-magnitude diagrams of ω Cen using JWST filters. (a) mF150W2 versus mF150W2 − mF322W2 CMD. (b) mF322W2 versus mF150W2 − mF322W2 CMD. In both panels, black dots represent stars passing the quality selections while grey dots represent stars not passing the selections. The shaded light red and dark red regions highlight areas of the CMDs affected by saturation. The light red region corresponds to saturated photometry, while the dark red region represents severe saturation, where the photometry is saturated even in the frame zero.

In the text

	Fig. 3. Same as Fig. 2, but only for the sources located in the overlap region between the JWST and HST FOV.
In the text

Fig. 4. Cluster membership selection based on PMs, for stars that passed the photometric quality selections and have measurable PMs. (a)–(b) VPD and mF322W2 versus mF150W2 − mF322W2 CMD, respectively. (c) mF322W2 magnitude versus the one-dimensional PM (μR). The red line separates cluster members from field stars. (d)–(e) VPD and mF322W2 versus mF150W2 − mF322W2 CMD for stars that passed the PM selection. (f)–(g) same diagrams for stars that did not pass the PM selection. In panels (c) through (g), sources that passed the PM selection are represented by blue dots, while those that did not pass are depicted as orange crosses. In panels (b), (e) and (g), the shaded light red and dark red regions highlight areas of the CMDs affected by mild and severe saturation, respectively.

In the text

Fig. 5. Crossing of the bMS and rMS stars in the JWST photometry-based CMD. (a) In the HST mF814W versus mF606W − mF814W CMD, we isolated a sample of bMS (black points within blue circles) and rMS (black points within red circles) stars; the two sequences are clearly distinguishable. (b) The same samples highlighted in the JWST mF150W2 versus mF150W2 − mF322W2 CMD, using the same colour code. The shaded light red and dark red regions highlight areas of the CMDs affected by mild and severe saturation, respectively.

In the text

Fig. 6. Comparison between the JWST NIR CMD of ω Cen presented in this study, and model isochrones. For clarity, we do not plot isochrones on the observed CMD directly. Instead, we compare the isochrones to the region of the colour-magnitude space occupied by the cluster, which is shaded in yellow in all panels. Panel (a) demonstrates how this region was chosen in relation to the observed CMD. Panel (b) over-plots two of the isochrones from Gerasimov et al. (2022) with both their original helium mass fraction (Y = 0.4, dashed lines) and with modified evolutionary models that incorporate a near-solar helium abundance (Y = 0.25, solid lines). Panels (c) and (d) show the Y = 0.25 NOM isochrone as well as its lower-MS variations with various chemical offsets. In the legend, ΔC, ΔO and ΔM refers to the offsets in [C/H], [O/H] and [M/H] from the chemical composition adopted for NOM in Gerasimov et al. (2022), where M only refers to metals heavier than oxygen. The shaded light red and dark red regions highlight areas of the CMDs affected by mild and severe saturation, respectively.

In the text

Fig. 7. Effect of offsets in the abundances of selected metals from the chemical composition of the NOM isochrone (panels (b)–(d) of Fig. 6) on the observed mF150W2 − mF322W2 colour of the star at Teff = 4000 K. The effects are shown at three different metallicity offsets, where M only includes metals heavier than oxygen. Positive values indicate that the star gets redder when the abundance of the element is increased. The 0 (no-effect) horizontal line is shown for reference.

In the text

Fig. 8. mF322W2 versus mF150W2 − mF322W2 CMDs of our selected sample of stars, focusing on the upper-MS (panel a) and the entire MS (panel b), with the two fiducials used to generate the ASs. The dashed grey rectangle in panel b highlights the CMD region displayed in panel a. The blue fiducial corresponds to the bMS, while the red fiducial corresponds to the rMS. In the upper part of the MS, the bMS occupies the blue side of the CMD, while the rMS is on the red side. Below the intersection point, the two sequences switch positions. The shaded light red and dark red regions highlight areas of the CMDs affected by mild and severe saturation, respectively.

In the text

Fig. 9. (a) Injected ASs for the two sequences (blue and red points) in the mF150W2 versus mF150W2 − mF322W2 CMD. Recovered ASs are represented in black. (b) Same as panel (a), but the black points now represent the real stars. (c)–(d) Completeness as a function of the mF150W2 and mF322W2 magnitude, respectively. The completeness for each sequence is represented using the same colour scheme as in panel (a). Shaded light red and dark red areas qualitatively mark magnitude ranges where MS stars are mildly or severely saturated in at least one filter.

In the text

	Fig. 10. Completeness as a function of mF150W2 (panel a) and mF322W2 (panel b) magnitudes for the two sequences in the overlapping region of the HST and JWST fields.
In the text

	Fig. 11. Same as Fig. 10 but including the PMs selection.
In the text

Fig. 12. Observed LFs and derived MF of ω Cen. Panels (a) and (b) show the observed LFs of ω Cen in F150W2 and F322W2 bands, respectively. The LFs are represented as histograms with 15 uniform bins. Over-plotted in both panels are the best-fitting LF models (taken to be the model with the smallest GoF in F150W2) from each of the five families listed in Table 1. The model curves include completeness effects. Panel (c) shows the MF of the cluster based on the inferred stellar masses from the observed F150W2 magnitudes and the Y = 0.25 NOM isochrone, shown in panel (b) of Fig. 6. Unlike panels (a) and (b), the inferred counts in panel (c) have been corrected for photometric completeness. Over-plotted are the best-fitting single-component and two-component power-law MFs. In the legends, pop stands for population, comp for component, and PL for power law. The light red areas qualitatively mark magnitude ranges where MS stars are mildly saturated in at least one filter. The LF and MF values in these regions should be interpreted with caution.

In the text

Fig. 13. Procedure for estimating the LFs of the sequences in ω Cen in the low-MS. Panel (a) shows the mF322W2 versus mF150W2 − mF322W2 CMD, focusing on the region where the separation between the components is most evident. The red and blue lines indicate the fiducials of the rMS and bMS, respectively, and are used to construct the verticalised CMD shown in panel (b). Panels (1,2c) through (1,2l) display the histogram distribution of ΔmF150W2 − mF322W2 for stars across eight magnitude intervals, as defined by the grey lines in panels (a) and (b). In panels (1c) to (1l), the red and blue regions corresponding to the two sequences are highlighted. These regions are used to estimate the LFs presented in Fig. 14, with values reported in Table 2. Panels (2c) to (2l) present the same histograms with the best-fitting multi-Gaussian models overlaid: three-Gaussian fits in panels (2c) to (2i) and a two-Gaussian fit in panel (2l). The individual Gaussian components are shown in red, green, and blue, while the combined fit is plotted in grey. The obtained LFs are shown in Fig. 15 and the values are reported in Table 3.

In the text

	Fig. 14. (a) LF derived using the regions defined in panels (1c) to (1l) of Fig. 13. (b) Population ratio of stars associated with the two sequences. In both panels, the red-shaded region indicates the saturation limit, while the grey-shaded region marks the threshold beyond which the three populations can no longer be distinguished.
In the text

Fig. 15. (a) LF obtained from the three-Gaussian fit shown in panels (2c) to (2l) of Fig. 13, highlighting the drop of the bMS in the final bin (mF150W2 ∼ 23.75). (b) Population ratio of stars associated with the three sequences. As in Fig. 14, both panels feature a red-shaded region indicating the saturation limit and a grey-shaded region marking the threshold beyond which the three populations can no longer be distinguished.

In the text

Fig. 16. (a) MFs of the three identified sequences, along with interpolating lines in the mass range 0.1 < M < 0.3 M⊙ (solid lines). The corresponding power-law slopes are indicated in the plot. For the gMS and bMS, we also show fits in the narrower mass range 0.2 < M < 0.3 M⊙ (dashed lines), with their respective slopes also reported. The black histogram represents the combined MF, as shown in panel (c) of Fig. 12. The solid black line and dashed yellow line correspond to the single-component and broken (two-component) power-law fits to the combined MF, respectively, with their slopes provided for reference. (b) Population ratio of stars associated with the three sequences as a function of stellar mass. The shaded light red and dark red regions highlight areas affected by mild and severe saturation, respectively.

In the text

Fig. A.1. Colour-magnitude diagrams of ω Cen using JWST filters. Panel (a) shows the mF150W2 versus mF150W2 − mF322W2 CMD, while panel (b) displays the mF322W2 versus mF150W2 − mF322W2 CMD, including all stars within the JWST FOV. In both panels, black dots represent stars that pass both the photometric quality and PM selection criteria, while grey dots indicate the remaining stars. The shaded light red and dark red regions highlight areas affected by saturation: the light red region corresponds to saturated photometry, while the dark red region marks severe saturation, where even frame zero is saturated.

In the text