© The Authors 2025

1 Introduction

Owing to photometric and spectroscopic observations conducted with various ground-based and space telescopes, recent studies have revealed that nearly all globular clusters (GCs) within the Local Group exhibit multiple stellar populations (MPs). The phenomenon of MPs is characterized by the fact that nearly all member stars within GCs exhibit abundance dispersions, involving light elements such as He, C, N, O, Na, Mg, and Al, and some even exhibit dispersions in heavy elements (see Gratton et al. 2004; Milone & Marino 2022, for an early and a more recent reviews). As expected, the phenomenon of MPs was first, and most extensively, reported in the GCs of the Milky Way (e.g., Smith 1987, for an early review). Subsequently, since last decade, it was discovered that star clusters with MPs also exist in the satellite galaxies, indicating that the MPs may be widespread among star clusters in different galaxies (Larsen et al. 2012; Milone et al. 2020, for examples in Fornax dSph and Magellanic Clouds).

The nature of the polluters remains debated, including asymptotic giant branch (AGB) stars (e.g., D’Ercole et al. 2008; Siess 2010; Doherty et al. 2014), fast-rotating massive stars (FRMS, Decressin et al. 2007; Krause et al. 2013), very massive stars (VMS, Vink 2018; Gieles et al. 2018), and interactive binaries (IBs, de Mink et al. 2009; Jiang et al. 2014; Nguyen & Sills 2024), or some of their combinations (e.g., Wang et al. 2020; Winter & Clarke 2023). For any scenario aimed at explaining elemental enrichment, given that helium is the most direct product of stellar nucleosynthesis, the polluters proposed by these models must also exhibit He-enrichment. Most importantly, these models would indicate different relationships between the dispersions of light elements and He abundance. Thus, studying He abundance distributions among different populations in GCs is crucial to understand the origins of multiple populations.

However, directly measuring He abundance of stars in GCs through spectroscopy is challenging. This is because He photospheric absorptions cannot be excited in the atmospheres of late-type stars, and the measurement of He abundance in stars with atmospheric temperatures exceeding 11 000 K is affected by the Grundahl jump effect (Grundahl et al. 1999), direct measurements of He abundance in GCs are very limited. For cool stars, typically of FGK type, only the chromospheric He I 10 830Å spectral line can be used for He abundance determination (e.g., Dupree & Avrett 2013). However, a recent study presented by Jian et al. (2024) has revealed that the strength of this absorption line strongly correlates with the Ca II log R′_HK index, indicating that the structure of the stellar chromosphere plays a significant role in the He I 10 830 Å line. This makes it inaccurate to characterize stellar He abundance solely using the He I 10 830 Å absorption. Consequently, in GCs, only horizontal branch (HB) stars within a specific temperature range can be employed for the determination of He abundance (e.g., Gratton et al. 2013; Marino et al. 2014). This makes it challenging to obtain a complete understanding of He abundance distribution among stars in GCs through spectroscopy. First, not all GCs have HB stars in this specific temperature range. Second, in GCs with MPs, the HB stars from the first and second generations differ significantly in color, making it difficult to measure the He abundance of HB stars for both populations. Because of this, in addition to direct spectroscopic observations, He variations among GCs are explored through photometric studies, such as extended or multiple HBs (e.g., Dalessandro et al. 2013; Milone et al. 2018), extended or dual red giant branch (RGB) bumps (Nataf et al. 2011), extended or multiple main sequences (MS, Norris 2004; Piotto et al. 2007; Cadelano et al. 2023) and their combinations (Bellini et al. 2013; Nardiello et al. 2015). Studies on helium abundance in clusters within satellite galaxies have only recently increased due to their great distance. Lagioia et al. (2019) reported evidence of He enhancement (δY ∼ 0.01) in the second population stars from four GCs in the Small Magellanic Cloud (SMC). Chantereau et al. (2019) found that four of the eight clusters in the Magellanic Clouds show indications of helium variation based on their HB morphologies. Our previous study has shown that the Large Magellanic Cloud (LMC) cluster NGC 2210 shows evidence of He variations of approximately δY =0.06–0.07 among its MS dwarfs (Li et al. 2023). In this study, we found that this research underestimated the impact of photometric errors in the cluster’s central region. We have corrected this value to at least (δY = 0.08), but this adjustment applies only to the cluster members of NGC 2210 outside the central area. Thanks to the Multi Unit Spectroscopic Explorer (MUSE) integral-field spectrograph based on the Very Large Telescope (VLT), it subsequently became feasible to perform spectral measurements of the MS stars within MC clusters, although measuring their He abundance can currently be made only for the early-type MS stars in young MC clusters. Carini et al. (2020) estimated the He abundance of ten bright MS stars in the SMC cluster NGC 330, and they found no evidence of variation in helium abundance among the individual stars. They subsequently analyzed another pair of double clusters in the LMC, NGC 1850A and NGC 1850B, and discovered significant differences in He abundance among the B-type main-sequence stars, with those belonging to NGC 1850B exhibiting a helium abundance exceeding that of NGC 1850A by more than (δY = 0.1, Carini et al. 2024).

Using the morphology of the HB to fit the He distribution of stellar populations is the most common method, as HB stars are very bright and their morphology often shows a wide range in color and magnitude, which helps reduce uncertainties caused by photometric errors. However, because HB is a very advanced phase of stellar evolution, models that rely on the analysis of HB morphology are usually subject to larger systematic errors. It is well known that the He variation is not the only parameter governing the HB morphology. Multiple factors may work together to influence the morphology of the HB, including the metallicity of stellar populations, age dispersion, and mass loss rates. Especially in GCs that contain MPs, HB stars from different stellar generations may undergo significantly different mass loss histories (Tailo et al. 2020). This further complicates the ability of the HB morphology to accurately represent the differences in their He abundances. Using the RGB bump to study the He abundance distribution in star clusters results in lower uncertainty due to the influence of stellar mass loss (e.g., Lagioia et al. 2018). However, this method is still affected by the metallicity of the stellar population. Moreover, even with the reasonable assumption that the metallicity of stars in the cluster is uniform, unresolved binary stars can substantially alter the luminosity function distribution of the RGB, thus obscuring the morphological features of the RGB bump (e.g., Li et al. 2023). Importantly, the RGB bump stars in a cluster typically make up a small proportion of the total population, making it nearly impossible to detect a clear RGB bump structure in GCs that lack a significant number of red giant stars.

Compared to the methods mentioned above, examining the morphology of the MS of GCs to find their He abundance distribution is more physically reliable. As helium is the second most abundant element in stars, variations in its abundance have a significant impact on stellar structure and evolution. First, the initial helium abundance in MS stars of a specific mass affects their effective temperature (T_eff) because the associated decrease in the radiative opacities of their outer layers makes them hotter as the He abundance (Y) increases. Second, an increase in the Y raises the molecular weight of the stellar gas, causing MS stars of a specified mass also to exhibit greater luminosity. Moreover, stars that are rich in helium evolve more quickly than normal stars, as they possess higher luminosities. The overall effect is that in the MS region below the MS turn-off (MSTO), He-enriched MS stars will appear significantly bluer compared to normal stars. As they approach the MSTO, the branch of He- enriched stars will gradually come closer to the normal MS, eventually connecting with the red end of the MSTO region and forming a dimmer subgiant branch (SGB). MS stars have not experienced significant mass loss or convective dredge-up, making their He abundance a more reliable indicator of the primordial He distribution of the GC. A series of studies have demonstrated that with the correct selection of optical bands, the multiple MSs seen in clusters can accurately reflect the internal He variations (e.g., Sbordone et al. 2011; Cassisi & Salaris 2013). Therefore, even though the luminosity of the MS is considerably less than that of the RGB bump and the HB, examining the He distribution in clusters by analyzing the morphology of the MS is a worthwhile method, especially when complemented by deep exposure observations, particularly from space telescopes.

To measure the He abundance of stellar populations through the morphology of the MS, the selected photometric bands must meet the following conditions: (1) The color should be sensitive to variations in helium abundance (i.e., dY/dcolor). (2) The color should not be sensitive to changes in other light elements due to the presence of MPs. (3) The relationship of dY/dcolor should account for different physical inputs (such as adopted mixing length αMLT , and/or atomic diffusion efficiency). A detailed theoretical study has shown that the F606W and F814W bands of the Hubble Space Telescope’s Advanced Camera Survey (HST/ACS) effectively satisfy the requirements for photometric passbands. The color distribution (F606W − F814W) of the MS obtained from these bands can reliably represent the He abundance distribution within stellar populations (Cassisi et al. 2017).

In this study, we analyzed the MS morphology of four GCs in the LMC, Hodge 11, NGC 1841, NGC 2210, and NGC 2257, with the aim of determining the He abundance distribution among their member stars. Our analysis relies on color-magnitude diagrams (CMDs) utilizing the F606W and F814W passbands of the HST/ACS. We observed a significant broadening of the MS in all four clusters, which, after minimizing the effects of field star contamination, photometric errors, unresolved binaries, differential extinction, and variations in metallicity, can only be explained by He abundance variation. We infer that there are likely significant variations in He abundance of at least δY = 0.08 within these clusters.

The following sections of this paper will include the subsequent topics: Section 2 will describe our methods for data reduction and our analysis, including photometry, cluster members determination, isochrones fitting, and other statistical analysis using ASs. Section 3 presents our main results, primarily focused on the He distribution of MS stars within the clusters, including the values of He enrichment (δY) of MS stars and the fraction of He-enriched stars. Section 4 contains a scientific discussion addressing these results along with our conclusions.

2 Data reduction and methods

2.1 Calculation of structural parameters

Our analyses are on four GCs within the LMC, Hodge 11, NGC 1841, NGC 2210, and NGC 2257. We first collected the spatial distribution of all stars with G ≤ 20 mag within 10 arcminutes of the centers¹ of these four clusters (∼ 144 pc at the distance of the LMC), utilizing the archival dataset from the Gaia DR3 (Gaia Collaboration 2023). We employed these catalogs to determine their structural parameters by fitting the radial profiles of their number density with a King profile (King 1962): $$\rho (r)=k{\left[\frac{1}{\sqrt{1+{\left(r/r_{\text{c}}\right)}^2}}-\frac{1}{\sqrt{1+{\left(r_{\text{t}}/r_{\text{c}}\right)}^2}}\right]}^2+b,$$(1)

in this equation, the core and tidal radii are referred to as r_c and r_t, respectively. The symbol b represents the number density of the background field, while k is a normalization coefficient. Additionally, ρ stands for the number density, and r indicates the distance from a star to the center of the cluster. The best-fitting King profiles are employed to calculate the half-number radii, referred to as r_h . For old GCs, the majority of stars are low-mass stars, which contribute most of the cluster’s total mass. Therefore, the half-number radius can generally be considered equivalent to the half-mass radius. The derived b are subsequently used for field decontamination. In Figure 1 we show the number density profiles of our clusters as well as the best fitting King profile. The derived structure parameters are summarized in Table 1.

We compared our results with those of Lanzoni et al. (2019) and found that, in terms of r_c and r_h (in their results, r_h means half-mass radius), our findings are consistent with their trends. For r_c, the size ranking of the four GCs is NGC 2210 < Hodge 11 < NGC 1841 ≈ NGC 2257. Similarly, for r_h, the size ranking again remains NGC 2210 < Hodge 11 < NGC 1841 ≈ NGC 2257. Within the margin of error, our r_c are systematically larger than those of Lanzoni et al. (2019). In contrast, for r_h, our results are consistently larger than theirs. Notably, our analysis suggests that Hodge 11 may possess the largest rt among the four clusters, whereas their fitting placed Hodge 11 third in the ranking. However, our fitting results for Hodge 11 exhibit the largest errors among the four clusters, which may be attributed to significant fluctuations in the stellar density surrounding Hodge 11 (see the top panel of Figure 1). Given that we utilized different datasets for our fitting and that our analyzed field of view (FoV) is larger than theirs (their fittings are based on the HST ACS/WFC data), the observed discrepancies are acceptable. Therefore, we conclude that our fitting is generally in agreement with the findings of Lanzoni et al. (2019).

Fig. 1 Stellar Number Density Profiles of Hodge 11, NGC 1841, NGC 2210, and NGC 2257. Best-fitting King profile shown by solid red curves. Solid blue line, dashed orange line, and dash-dotted green line indicate core radius, half-number radius, and tidal radius, respectively.

2.2 Photometry

Our analysis mainly used data obtained from HST ACS/WFC observations because it has a deeper exposure than the Gaia DR3. The program ID is GO-14164 (PI: Sarajedini). To improve the accuracy of removing field star contamination, we also used the observational data from the Ultraviolet and Visible Channel of the Wide Field Camera 3 (UVIS/WFC3) on the HST for the four clusters in a similar F814W passband. This program ID is GO-16748 (PI: Niederhofer). The average timeline between the two observations for our four clusters is ∼6 yr. We will subsequently employ both datasets to obtain the relative proper motions (PMs) of the stars within their overlapping FoV, thereby facilitating a more precise field star decontamination. The basic information about the dataset we used is summarized in Table A.1.

For each cluster, we use the specific HST photometric package DOLPHOT2.0 (Dolphin. 2011a,b, 2013) to conduct pointspread function (PSF) based photometry on their charge-transferefficiency-corrected frames (the images with root names ending in ‘_flc’ and ‘_drc’). The photometric procedure we employ involves the following steps. First, we identify and mask bad pixels in the images. This step is automatically carried out using the acsmask command for ACS/WFC data and the wfc3mask command for UVIS/WFC3 data. Next, we separate the chargetransfer-efficiency-corrected frames into two chips using the splitgroups command. We then compute the corresponding sky images for all frames in order to estimate the background levels at different locations, utilizing the calcsky command. In this step, we set the parameters based on the tutorial’s recommendations, choosing step=4 in the sky-fitting settings. At this point, the program will fit the sky within the photometry aperture as a four-parameter PSF fit, which is the most accurate sky calculation provided by the photometric package. After completing the previous steps, we begin photometry using the dolphot command. At this stage, we will choose the drizzled image (‘*_drc.fits’, per filter) as the reference frame. Since we will attempt to reduce field star contamination by using the changes in star positions, PM, from two observations, we need to eliminate errors caused by factors such as CCD distortion and alignment as much as possible. We use the WCS header information for alignment, setting UseWCS=2 to estimate a full distortion solution.

The original star catalog we obtained contains a large amount of information, including a total catalog formed by merging all frames and individual catalogs from each separate frame. Since we are analyzing MS stars, short exposure times can lead to high photometric errors, making it impossible to estimate the distribution of He abundance within the cluster. Therefore, we decide to choose the combined catalog from all frames as our analysis sample. We first remove false detections identified by DOLPHOT; these detections will have a magnitude of 99.999 mag in at least one passband. We then carried out a quality filter on the valid detections, choosing only the stars clearly labeled as good stars (Object types=1). We removed sources with sharpness values less than −0.3 and greater than 0.3. According to the DOLPHOT manual, these sources are more likely to be cosmic rays (>0.3) or extended sources (such as blends or galaxies, <−0.3). Roundness describes the ellipticity of a source’s light profile. A perfectly circular source has a Roundness value of 0. In this study, we exclusively select sources with a Roundness value lower than 0.5, as sources exhibiting higher roundness are unlikely to be perfect stars (e.g., those blended with diffractive spikes). Crowding refers to the level of contamination caused by nearby sources. The Crowding value, expressed in magnitudes, indicates how much brighter the star would have appeared if the adjacent stars had not been fitted simultaneously. We applied a fairly relaxed limitation on crowding, excluding only sources with a crowding value greater than 0.5 mag. We discovered that most of the MS stars we are interested in had significant crowding, as their intrinsic brightness was similar to the background level. Finally, we also removed sources with an SNR < 5, as these sources had too large errors that could significantly affect the analysis. We will then analyze artificial stars (ASs) that went through the same selection criteria to statistically compare whether the morphology of the MS matches the predictions of the simple-stellar population (SSP) models.

2.3 Estimation of stellar proper motions

The header files of the images used for photometry contain World Coordinate System (WCS) keywords, which define the relationship between pixel coordinates in the image and celestial coordinates. Specifically, these files include eight parameters, among which CRPIX1 and CRPIX2 represent the pixel coordinates of the reference point to which the projection and rotation pertain. CRVAL1 and CRVAL2 specify the central coordinates in terms of right ascension (RA, or α) and declination (Dec, or δ), or longitude and latitude, expressed in decimal degrees. The FITS WCS standard employs a rotation matrix consisting of CD1_1, CD1_2, CD2_1, and CD2_2 to indicate both rotation and scale, facilitating more intuitive computation in cases where the axes are skewed. In particular, the conversion of the X and Y coordinates on the CCD to RA and Dec can be done using the following relationships: $$\begin{array}{l}\delta =\text{ CRVAL}2+\left[\text{ CD}2\_1(X-\text{ CRPIX1 })+{\text{CD2}}_{\_}2(\text{Y- CRPIX 2})\right]\\ \alpha =\text{ CRVAL}1+\frac{\left[\text{ CD1}\_1(X-\text{ CRPIX}1)+\text{ CD}1\_2(\text{Y}-\text{CRPIX}2)\right]}{\cos \delta }.\end{array}$$

To more accurately identify the member stars of the star cluster, we used a “fool man’s method” to assess the PMs of each star using datasets with a time interval of about 6 years. For each star observed by ACS/WFC around ∼2016 (first observation), we looked for its matching observation in the UVIS/WFC3 data from ∼2022 (second observation) by pairing it with the source that was nearest in distance during the second observation. For each star, we recorded the distance to the nearest star identified in the ∼2022 observation. We found that the distribution of the nearest distances showed three clear different peaks. The first peak represents stars that were not covered by the FoV of the second observation due to incomplete overlap between the two observational fields. In this situation, our method randomly selects a star from the second observation as its nearest counterpart, leading to very large nearest distances for these stars, which we then excluded first. The second peak mainly results from the different exposure depths of the two observations, which causes some stars to be missed in the second observation. In this situation, our method still randomly picks the nearest star from the second observation as its match; however, since both stars are now in the same FoV, the nearest distance is much smaller than in the first peak. Nevertheless, we can easily remove these stars by looking at their two-dimensional distribution in (∆α) and (∆δ). The third peak corresponds to nearest distances that are much smaller than those of the first two peaks, representing the same stars detected in both observations. From the distribution, this nearest distance is very close to zero, as the movement of the stars over six years is not enough to create noticeable changes in position. However, we still see that this peak has an extension towards larger distances, which we think comes from the influence of field stars. To more accurately distinguish between field stars and member stars, we further analyzed the distribution of these stars in a two-dimensional PM space (∆α cos δ, ∆α)). We found that their overall distribution remains relatively continuous, and it does not exhibit an isotropic uniform distribution. This indicates that relying solely on a 6 yr time baseline results in low accuracy for measuring PM in terms of direction. We then created iso-density contour maps of these stars based on their distribution in the two-dimensional PM diagram. Using the previously calculated field star contamination rate, we chose specific contour lines as boundaries. The stars inside these contour lines are classified as member stars of the star cluster, while those outside are considered field stars.

The PM distribution of these stars and the contour lines we selected are displayed in the subplots of Figure 2, and Figures B.1–B.3. These four figures show the CMDs of all observed stars in four clusters (left panel), the CMDs of the selected cluster member stars (middle panel), and the CMDs of the identified field stars (right panel). In the upper right corner of each diagram, we provide two subplots: the PM distribution of the corresponding stars (top) and their spatial distribution (bottom). We found that, based on the contamination rate we calculated and the relative PM distribution from two observations, we could effectively remove most field stars. Specifically, the CMDs of the removed field stars are quite scattered, clearly deviating from the isochrones that describe the star cluster. Additionally, the spatial distribution of the removed field stars is mostly uniform, showing no clear clustering, indicating that we did not overestimate field star contamination² (Except for NGC 2210, we found that no matter how we adjusted the criteria for removing field stars, the eliminated field stars still showed some clustering. In the future, with updated observational data (allowing for PM calculations based on a longer age baseline), we may be able to improve this situation).

Fig. 2 CMD of Hodge 11. Left: CMD shows all stars observed in the field of cluster Hodge 11. Middle: CMD displays stars classified as cluster members. Right: CMD features stars identified as field stars. In each panel, inset plots illustrate differences in coordinates between two observations taken approximately 7 years apart (top) and the spatial distribution of the corresponding stars.

2.4 Artificial stars

DOLPHOT allows for the creation of ASs on the original images, based on the PSFs used in the initial photometry. Before this, we need to set the magnitudes of the ASs. To make the ASs as realistic as possible, we first fit isochrones to the observed star cluster. The isochrones we used are derived from the Princeton- Goddard-PUC (PGPUC) stellar evolutionary code³ (Valcarce et al. 2012, 2013). We used this model because it is specifically made for studying low-mass dwarfs, which perfectly fits our research goals. The inputs for the PGPUC model comprise the values of age, He mass abundance (Y), global metallicity (Z), solar-scaled abundance of α-element ([α/Fe]), and mass-loss rate (η). Among the parameters mentioned, age, Y, Z, adopted distance modulus, (m − M)₀, and extinction, E(B − V), have the most substantial impact on our fitting, while [α/Fe] and η exert minimal influence. For all star clusters, we have consistent inputs of [α/Fe]=+0.3 (Wagner-Kaiser et al. 2017) and η=0.2 (default values). The best-fitting isochrones are determined through visual inspection. As a result, our estimation of the errors in the fitting parameters is rather informal, simply matching the grid density used during our fitting. Furthermore, Y =0.25 is always kept fixed in our fitting, as it serves as a reference for creating MPs with different Y values. The fitting parameters are summarized in Table 2. Our fitting results are consistent with those of Wagner-Kaiser et al. (2017) within the margin of error.

Clearly, using visual inspection for isochrone fitting introduces considerable subjective bias. To create synthetic stellar populations that closely match real observations, we use the MS ridgeline (MSRL) as a reference for generating ASs. Similar to Li et al. (2023), we determined the MSRL below the TO region using the robust regression method based on the Gaussian process and iterative trimming created by Li et al. (2020, 2021).

The PGPUC isochrones in our analysis mainly aim to show how the colors of stellar populations change with varying Y, i.e., ∆Y/∆(F606W − F814W). We subsequently calculated 11 isochrones with varying He enrichments of ∆Y = 0.01 to 0.12 (Y = 0.26 to 0.37, with an increment of 0.01). Thus, for each star cluster, we calculated 13 synthetic stellar populations with different values of Y (including the fitted value of Y = 0.25). The color deviations for each isochrone relative to the standard isochrone (Y = 0.25), represented as ∆(F606W − F814W), are then incorporated into the MSRL.

We then used DOLPHOT to generate ASs. The magnitudes of the ASs are derived from the MSRL, while their luminosity function is derived from the standard Kroupa mass function. Each population consists of 2 million ASs. Therefore, for each cluster, we generated 2.6 × 10⁷ ASs with Y values from 0.25 to 0.37. Since ASs are created using the same PSF as real stars, and we have filtered them as we did for real stars, they can, in principle, simulate an observation with similar photometric uncertainties (including noise from cosmic rays, hot pixels, crowding effects, and other artifacts) as real stars. For each synthetic population, we added additional color noise to simulate the residuals after the removal of differential reddening. The differences between the recovered and input magnitudes of ASs are referred to as the photometric uncertainties⁴.

Based on the photometric errors derived from ASs, we decided to analyze MS stars within the magnitude range of 22.5 mag < F606W < 24.5 mag. MS stars brighter than this magnitude range are too close to the MSTO region, where increased helium abundance accelerates stellar evolution, which offsets the effects of luminosity and temperature changes caused by He enrichment. Below this magnitude range, MS stars have too low SNRs, leading to the spread of the MS being dominated by photometric uncertainties, which would mask any possible signals of He variation. For each cluster, ASs indicate an average completeness of ≥65% for MS stars within this magnitude range.

Through ASs, we also found that, within ∼38 arcsec (about 9 pc) from the center of each cluster, the photometric uncertainties increase sharply due to crowding. Therefore, we decided to analyze only the sample of MS stars located between 10 pc and 30 pc. Beyond 30 pc, the proximity of stars to the edge of the FoV also leads to a slight increase in both photometric and astrometric uncertainties, see Figure C.1.

2.5 Differential extinction

The dust distribution in the ACS/WFC FoV (202″× 202″ , corresponding to 48.5 pc × 48.5 pc at the distance of the LMC) may be inhomogeneous. We have used the method outlined in Milone et al. (2012) to minimize the effect caused by the dust inhomogeneity – the differential reddening. Briefly, to quantify the amount of differential reddening suffered by each star across the FoV, we chose RGB and SGB stars instead of MS stars as the reference population to avoid underestimating the helium spread. We derived the fiducial line of the reference population based on a 2D probability density function in the F606W versus F606W-F814W CMD. To simplify the correction process, we rotated the CMD counterclockwise to define a new reference frame in which the abscissa (X′) and the ordinate (Y′) are parallel and orthogonal to the reddening direction, respectively. For each star, we calculated the average systematic offset ($\overline{\text{Δ}{X}^{\prime }}$) from the fiducial line, based on the nearest 50 reference stars, excluding the target star from the computation. To correct the photometry for the effects of differential reddening, we subtracted the derived $\overline{\text{Δ}{X}^{\prime }}$ from the X′ value of each star. After applying this correction, the de-reddened CMD was used to refine the selection of reference stars and update the fiducial line. This process was iterated until convergence, typically after about three iterations. The corrected X′ and Y′ were transformed back into F606W and F814W magnitudes, using the relative absorption coefficients from the extinction curves of Cardelli et al. (1989) and O’Donnell (1994).

We constructed a series of extinction distributions utilizing ASs (as detailed in the following subsection) and then calculated the extinction of the ASs using the same method employed for real stars. By comparing the calculated extinction to the input values, we estimate the extinction residuals associated with our approach. Specifically, we randomly selected ten points within the FoV of the ASs, assuming their extinctions were randomly drawn from a distribution function with varying dispersion. We then performed a two-dimensional polynomial fitting of extinction based on these ten points across the entire FoV. The extinction for other ASs in the field was obtained by interpolating this fit, resulting in a specific differential extinction distribution. We applied the input extinction to correct the input magnitudes of each AS by using the same method applied to the real observation, and obtained the photometric magnitudes of the ASs. We used the above statistical methods to calculate the differential extinction for each AS, and the calculated differential extinction was used to correct the photometric magnitudes (not input magnitudes) of the ASs. We also performed photometry on the ASs that did not include the differential extinction. By comparing the ASs that included differential extinction with those that did not, we obtained the color dispersion difference in their MS, δ(F606W − F814W), which we conclude describes the extinction residuals after de-reddening.

In Figure D.1, we present an example to demonstrate the accuracy of our de-extinction method. The upper left panel shows a spatial distribution of an AS population, with color representing their input extinction. The upper right shows the same distribution but with color representing their calculated extinction. We found that the calculated extinction distribution matches the input distribution well. The bottom panel shows the correlation between the input differential extinction and the color residuals after correcting for differential extinction. We found that, except for one outlier, there is a positive correlation between the input differential extinction and the color residuals. We use a linear function to fit this positive correlation and make corrections to the subsequent analyses of the ASs.

2.6 Quantify helium spread

Our subsequent analysis is essentially the same as in Li et al. (2023). For both the observed and simulated MSs, we calculated the color differences, ∆(F606W − F814W), for each star relative to their MSRL. Here the simulated MSs can be composed of either a single population (i.e., an SSP with a fixed Y value) or MPs with different Y values. From each artificial stellar population, we select a subsample that matches both the luminosity function and the total number of stars of the real observations, serving as a representative sample. A synthetic MPs is then constructed as a composition of these synthetic populations. To minimize the influence of unresolved binaries, we selected only those MS stars located on the blue side of the MSRL for our analysis. For the sake of clarity, we “folded” the observed and simulated MS stars on the blue side of the MSRL, assuming that there exists a corresponding mirror star at a symmetrical position relative to the MSRL. We divided the observed and simulated MS stars into 76 bins based on their ∆(F606W − F814W) values (where the number of bins is roughly equal to the square root of the number of observed stars) and counted the corresponding number of stars in each bin for both the observed and simulated data. We used the χ² minimization method to evaluate the similarity between the observed and simulated distributions: $${\chi }^2={\displaystyle \sum _i\frac{\left(N_{\text{i}}^{\text{obs}}-N_{\text{i}}^{\text{mod}}\right)}{N_{\text{i}}^{\text{obs}}}}$$(2)

where $N_{\text{i}}^{\text{obs}}$ and $N_{\text{i}}^{\text{mod}}$ are the number of observed and simulated stars with their relative colors in the i-th bin.

Our model contains two parameters: one is δY, defined as the difference between the maximum Y in the synthetic stellar populations and Y = 0.25 (the He abundance of the primordial stellar population). The other parameter, f_2P , represents the ratio of stars with Y > 0.25 to the total number of stars in the synthetic population. Similar to Li et al. (2023), we assume two extreme cases for the Y distribution: one where δY is uniformly distributed, meaning that the number of stars with He enhancement is the same across different Y values (for Y > 0.25). For example, if we set δY = 0.05, the model stellar population would contain six populations with Y values of 0.25, 0.26, 0.27, 0.28, 0.29, and 0.30. If we assume f_2P = 0.5, then the population with Y = 0.25 would represent 50% of the total number of stars, while the five He-enriched populations with Y values of 0.26, 0.27, 0.28, 0.29, and 0.30 would each makeup 10% of the total stellar count. The other case assumes a highly discrete distribution of δY, where the model MPs consist of only two populations, corresponding to Y = 0.25 and Y = Y_MAX. For example, in the case of δY = 0.05, the model would include only the populations with Y = 0.25 and Y = 0.30. If f_2P = 0.5, this means that both groups represent 50% of the total number of stars. The two scenarios described above represent two extreme distribution cases for a given f_2P . Therefore, the actual distribution of δY is likely to lie somewhere between these two cases.

Fig. 3 χ2 contour maps as a function of internal He spread, δY, and fraction of He-enriched stars, f2P . χ2 values are coded by colors. Best-fitting points (minimum χ2) are represented by red stars. Black solid lines indicate 1σ uncertainties, calculated from χ2 variations. The adopted model assumes a uniform He distribution. Panels a, b, c, and d correspond to clusters Hodge 11, NGC 1841, NGC 2210, and NGC 2257, respectively.

Fig. 4 Same as in Fig. 3, with the adopted model assumes a bimodal He distribution.

3 Main results

Our main results are shown in Figures 3 and 4. The two contour figures show the distribution of χ² obtained from comparing observations with models, with the value indicated by color. The x-axis represents the maximum He dispersion, δY, while the y- axis shows the ratio of stars with He enrichment to the total population, f_2P . The point with the minimum χ² is marked on the plot as a red pentagram. We fitted a third-order polynomial to the distribution of χ²(δY, f_2P) and calculated the residuals δχ² . The contour line corresponding to the minimum χ² + δχ² is marked with a solid black line, indicating the fluctuation level of χ² due to the random sampling effects of the model.

Figure 3 presents the results based on the assumption that δY follows a uniform distribution. We found that the derived best-fitting values for δY and f_2P are substantially high for all clusters. With the exception of NGC 1841, our fitting indicates that both Hodge 11, NGC 2210, and NGC 2257 exhibit an internal helium dispersion of at least (δY = 0.12), with approximately 80% of their stellar populations being He-enriched ( f_2P ∼ 0.8. In comparison to the other three clusters, NGC 1841 exhibits a slightly lower internal He dispersion and fraction of He-enriched stars, with δY =0.11 and f_2P=0.6. This implies that assuming a uniform distribution of helium, our fitting suggests that the stellar populations of all clusters are predominantly dominated by second-generation stars enriched in helium.

If we assume a bimodal distribution of helium, our approach would yield a relatively moderate He-enrichment among the stellar populations, which is present in figure 4. In this case, we find that the best-fitting (δY) values for the four clusters are between 0.08 and 0.10, and the fraction of He-enriched stars is about 40% to 50%. This means that the numbers of first-generation and second-generation stars are quite similar in these clusters. We present the details of our χ² fitting results in Table 3.

To guarantee the reasonableness of our results, we carefully selected several sets of best-fitting models and compared them with the observations through visual inspection. The results are illustrated in Figure 5. We discovered that the observed widths of their MSs are indeed all significantly larger than what is predicted for an SSP. This broadening cannot be explained by photometric errors, differential extinction (including residuals after de-reddening), or unresolved binaries. Visual inspection indicates that the observed MS aligns more closely with the morphology of MPs containing helium dispersion. If we assume a uniform distribution of helium abundance, larger values of δY and f_2P are required to reproduce the observed broadening of the MS. Conversely, if the distribution of helium abundance is bimodal, the observations can be reproduced with smaller values of δY and f_2P .

We employed the Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011) to investigate whether our findings are consistent with the HB morphology. Given the extensive computational time necessary for stellar evolution calculations using MESA, we performed a qualitative analysis to determine the δY required for a morphological fit to the observations. Therefore, the stellar populations we examined all comprised a similar number of stars, with the Y distribution ranging from 0.25 to 0.37, with increments of 0.01. The initial metallicity, Z_init, that we input to MESA is consistent with the parameters used for isochrone fitting, as outlined in Table 2. As described in Section 1, the morphology of the HB is significantly influenced by the stellar mass loss rate, η. This mass loss rate, in turn, depends on the metallicity of the cluster and the generation of the stellar population to which the stars belong (Tailo et al. 2020); thus, it remains highly uncertain in our simulations. The four-star clusters we examined have metallicities ranging from Z =0.0002 to 0.0005, corresponding to [Fe/H] ∼ −1.9 to −1.5. According to Tailo et al. (2020), first-generation stars within this metallicity range exhibit mass loss rates, η ∼ 0.25 to 0.35. Consequently, for our calculations, we simplistically adopted η=0.3, solely to qualitatively analyze the differences in He abundance distributions derived from the HB compared to those obtained from the MS.

Our results are shown in Figure 6, and we find that the He abundance distributions obtained from the HB of the four-star clusters do not always match those obtained from analyzing the MS. By visually comparing the observed HB morphology with the predictions from our models, we discovered that the two clusters, Hodge 11 and NGC 1841, may exhibit helium spreads of approximately δY ∼ 0.12 and 0.09, respectively. This is consistent with the He distribution obtained from analysis of their MS morphology (δY ∼ 0.10–0.12 and δY ∼ 0.08–0.11, respectively). However, the distributions of δY obtained from fitting the HB morphology of the other two clusters, NGC 2210 and NGC 2257, are very different from the results obtained from fitting the MS. The fitting results of the HB indicate that the helium abundance variation in NGC 2210 and NGC 2257 is around δY ∼ 0.03–0.05, which is only one-third to half of the results obtained from fitting the MS morphology. It is clear from Figure 6 that the HB of NGC 2210 and NGC 2257 is noticeably shorter than that of Hodge 11 and NGC 1841. Moreover, Figure 5 shows that the MS width of NGC 2210 and NGC 2257 is similar to Hodge 11 and is even a bit wider than that of NGC 1841’s MS.

4 Discussions and conclusions

We first discuss the reliability of our results. Many factors can contribute to the broadening of the MS. Besides the helium abundance dispersion, observational factors include photometric errors and PSF residuals. Physical factors include unresolved binaries, line-of-sight blending caused by crowding, differential extinction, dispersion of light elements (usually C, N, and O), metallicity dispersion, and field contamination. We discuss them separately.

DOLPHOT automatically provides magnitude uncertainties; however, these uncertainties are often inaccurate as they only account for the magnitude uncertainty arising from the counting uncertainty. For fainter sources, count uncertainty is the dominant source of error, while for the majority of brighter sources, it typically underestimates the true observational errors. Using ASs to simulate stellar populations can effectively avoid underestimating the errors mentioned above. This is because ASs are created based on the PSF fitted to observations, and their magnitudes are measured using photometry that matches the real observations. Therefore, in principle, ASs can account for errors from instrument effects such as bad pixels or saturation. Additionally, we add ASs in batches to the original frame for photometry, which allows us to consider errors from crowding, cosmic rays, and variation of background. In fact, by comparing the measured magnitudes of ASs to the input magnitudes, we find that the discrepancies are indeed greater than the theoretical errors provided by DOLPHOT.

PSF residuals not only depend on the PSF used, but often differ across different passbands, leading to effects similar to differential extinction and causing color spread to be greater than theoretically expected. However, since the direction of broadening caused by PSF residuals is not perpendicular to that of differential reddening, correcting for differential reddening automatically accounts for zero-point variations caused by PSF residuals. Therefore, we did not further iterate in the color direction to correct for the PSF residuals. We compared the results of correcting only for differential reddening and only for PSF residuals, and found that the color dispersion we obtained was nearly identical in both amplitude and spatial distribution.

The observational effects caused by unresolved binaries and blending are essentially the same (if both sources of blending are stars). ASs can statistically simulate the blending effect, but the impact of unresolved binaries cannot be eliminated. However, since both blending and unresolved binaries only shift photometric results towards the red end of the MS, analyzing only the blue half of the MS can help minimize the influence of binaries and blending.

If MPs are present, different stellar populations may have different light elements. Apart from He, the greatest variations are often seen in C, N, O, and Na, and changes in these elements can significantly affect photometry. According to Cassisi et al. (2017), C, N, O, and Na can produce significant bolometric corrections in the blue and UV passbands, while their impact in the optical bands is moderate, especially having an almost negligible impact on the F606W-F814W color band.

Through the fitting of isochrones, we can rule out the presence of metallicity spread, for reasons similar to those presented in Li et al. (2023). Changes in [Fe/H] would not only affect the dwarf stars but would also impact the colors of giants. Nevertheless, we do not observe any significant broadening of the RGB. Furthermore, we confirm that variations in [α/Fe] have nearly no impact on photometry.

We argue that the largest uncertainty in the broadening of the MS comes from field star contamination. Because the Magellanic Clouds are very far away, past attempts to decontaminate field stars have mostly used statistical methods. Specifically, this involves selecting a reference field in the outer areas of the FoV, far from the cluster center, or performing a parallel observation that is far away from the cluster involving the same passbands. It is assumed that the stars in these outer regions are mostly field stars. Then, a comparison of the CMDs between the outskirt and center regions of the cluster is made, followed by a statistical random removal of a certain number of stars from different parts of the diagrams (e.g., Li et al. 2016). Although this method can in principle statistically estimate the shape of the cluster’s CMD, it cannot determine whether specific stars are members of the cluster, which often leads to controversy (e.g., Cabrera-Ziri et al. 2016). Using stellar PM to determine whether a star belongs to the cluster is the only way to avoid such controversy (e.g., Wang et al. 2024).

We have confirmed that the use of previous statistical methods for field star decontamination does not significantly change the result that the MSs of the four clusters are significantly broadened (see, Li et al. 2023). However, the seven-year observation interval may only allow for a limited distinction between cluster member stars and field stars. Therefore, we expect that the so-called ‘cluster sample’ may still contain some field stars. Our fitting of the radial number density in the field of these clusters indicates that the contamination rate from field stars does not dominate. They only account for only 5% to 25% of the entire population. Therefore, even if field stars are not removed, we should expect at least three quarters of the dwarfs would distribute within a narrow MS, a typical SSP signal that cannot be obscured by a small number of field stars. However, the MS we observed is nearly uniformly distributed over a wide color range, suggesting that field star contamination alone cannot account for the observed MS morphology.

We acknowledge that some of our earlier conclusions about NGC 2210 are no longer reliable because we did not recognize the significant photometric errors in the central region of the cluster. Although we still stand by our main conclusion that NGC 2210 shows a significant broadening of the MS, our claim that its central region may possess a higher proportion of 2P stars is inaccurate (Li et al. 2023, their Fig. 9). Our conclusion was based on the observation of a greater number of MS stars with significant color dispersion in the center of the cluster, which we interpreted as indicating higher He abundance. However, this may simply result from larger photometric errors in the stars of the central region. We apologize for this error. Although in this work, we have excluded the very central region with larger photometric errors and only analyzed the relative outskirts of the cluster (>10 pc to the cluster center), where the photometric errors remain moderate and stable. However, this does not mean that crowding has no effect in this region; it only indicates that the contribution of crowding has become stable. The δY obtained by our method might only represent an upper limit. Similar to Li et al. (2023), we present our results using Galactic GCs (Milone et al. 2017a), two younger LMC clusters (Li et al. 2021; Ji et al. 2022), and SMC clusters (Chantereau et al. 2019; Lagioia et al. 2019) on the cluster mass − δY diagram. From Figure 7, it is evident that the four clusters we measured appear as outliers in the mass − δY diagram, independent of whether we use their present masses (top panel) or initial masses (bottom panel, see Baumgardt et al. (2019) for the calculation of initial masses). Next-generation instruments with higher spatial resolution may further refine our results.

The differences in results obtained from fitting the morphology of the HB and the MS clearly highlight the importance of mass loss. For example, the mass loss rate varies among different clusters and also between 1P and 2P population stars, even within the same cluster (Tailo et al. 2020). If a star is in a close binary system, the companion’s tidal forces will also affect its mass loss rate. A rapidly rotating star (for example, due to tidal locking that enhances the spin of a red giant) may have an extended red- giant lifetime, leading to greater mass loss, which results in its position in a bluer end of the HB. Overall, the complex evolutionary pathways of stars can result in different luminosities and colors during their HB phase, suggesting that the characteristics of the HB should not be the only criterion for evaluating a star’s He abundance. We also emphasize that Hodge 11 may be a particularly interesting LMC GC similar to the Galactic GC, NGC 2808, as its extended HB and broadened MS are consistent well with models exhibiting significant He enrichment (δY > 0.10).

If the estimated He abundance can represent the average level of the entire clusters, then the He enrichment in these four LMC clusters is comparable to that of the most massive GC in the Galaxy (>10⁶ M_⊙). This is surprising, as the current masses of these clusters are between 1–2.5 × 10⁵ M_⊙ (McLaughlin & van der Marel 2005). The f_2P for these four clusters are quite similar, therefore we did not observe a clear correlation between f_2P and their present-day masses. This is similar to the results found by Milone et al. (2020) in intermediate-age clusters of the Magellanic Clouds (illustrated in their Figure 7, top left). However, this outcome might also stem from the small size of our cluster sample, which could introduce variability. The derived f_2P for our four clusters ranges from ∼40–90%, depending on the chosen He distribution models. If present, these values would be comparable to those observed in Galactic GCs and are higher than those observed in intermediate-age clusters in the Magellanic Clouds, as reported by Milone et al. (2020). These four clusters are dominated by 2P stars, similar to the fact that extreme He-enriched clusters in the Milky Way are primarily occupied by 2P stars. This also differs from the old GCs in the SMC, where first-population stars are more dominant in number (e.g., Dalessandro et al. 2016; Lagioia et al. 2019). This may suggests that the MPs in these four LMC clusters may have origins similar to those of Milky Way GCs Milone et al. (2018), which is hard to explain because LMC clusters likely had weaker interactions with their host galaxy, keeping most of their original mass. This makes it unlikely for them to form such high fractions of second-population stars (Parmentier 2024). Again, given the potential effects of crowding, it remains challenging to assess the extent to which this result requires further refinement, however.

Assuming that our results accurately represent the helium spread within the cluster, they may provide insights into various MPs models. For the bimodal He distribution, our results indicate that the δY for the four clusters ranges from 0.08 to 0.10. Current mainstream models can generally produce such a significant He enrichment; however, VMS models tend to yield smaller helium enrichments (with δY < 0.01). The observed δY may suggest that the wind material from VMS originates from their later evolutionary stages (Gieles et al. 2018, their Section 4.3.2). However, if the actual He distribution is uniform, our results indicate that at least three clusters (Hodge 11, NGC 2210, NGC 2257) would have δY of at least 0.12. This approaches the maximum He enrichment limit that AGB models can produce (0.36–0.38, Siess 2010; Doherty et al. 2014). This suggests that AGB stars are unlikely to be the sole source of the observed He-enriched stars, or at least not the only contributors. FRMS models can generate significant He enrichment, reaching up to Y =0.8 (Decressin et al. 2007), while massive IBs can, in principle, produce a wide range of He enrichments by adjusting binary parameters (de Mink et al. 2009). If there are extreme He-enriched stars with δY > 0.12, we may need to consider the contributions from FRMS or IBs.

In conclusion, we find that the LMC GCs Hodge 11, NGC 1841, NGC 2210, and NGC 2257 have MS that is significantly wider than expected from the SSP model, which cannot be explained by photometric uncertainties, field contaminations, unresolved binaries, differential reddening and metallicity spread. We suggest that they may all have significant He enrichment among their MS stars. If present, this enrichment is similar to what is observed in Galactic GCs that exhibit extreme variations in He abundance. However, owing to the great distance of the LMC, our results are likely to represent an upper limit.

Fig. 5 Comparison of MS observations and synthetic MS. From left to right are comparisons between the observed MS, and the synthetic MS of the SSP, of the best-fitting MPs with a uniform He distribution, and of the best-fitting MPs with a bimodal He distribution. From top to bottom are clusters Hodge 11, NGC 1841, NGC 2210 and NGC 2257.

Fig. 6 Comparison of observed and simulated HB. From left to right are the observed and simulated HBs for clusters Hodge 11, NGC 1841, NGC 2210 and NGC 2257. color bars indicate the helium mass fraction, Y, of the simulated HB stars.

Fig. 7 Internal helium spread and cluster mass relationships. Top: the internal helium spreads and the clusters present-day masses diagram, δY − log( M/M⊙). Red pentagrams are LMC clusters. Small grey dots are Milky Way GCs. Dark/light grey circles are SMC clusters. Bottom: the correlation between the internal helium spreads and the clusters initial masses, for Milky Way GCs, LMC and SMC clusters. The arrow indicates that we consider this value to represent only an upper limit. The arrow indicates that we consider this value to represent only an upper limit.

Acknowledgements

This project is supported by the National Natural Science Foundation of China (NSFC grant Nos. 12033013 and 12073090), the National Key R&D Program of China (2020YFC2201400) and the China Manned Space Project with Nos. CMS-CSST-2021-A08, CMS-CSST-2021-B03.

Appendix A Datasets

This section provides detailed information about the HST datasets used in this work. Table A.1 includes the observational data we used, including the bands, exposure time, observation time, and Program ID.

Appendix B Field decontaminated color-magnitude diagrams

This section presents the field decontaminated CMDs of three other GCs: NGC 1841, NGC 2210, and NGC 2257. The decontaminated CMD of Hodge 11 is provided as an example in Figure 2, and relevant descriptions of these CMDs can be found in Section 2.3.

Fig. B.1 Same as in Fig. 2, but for the cluster NGC 1841.

Fig. B.2 Same as in Fig. 2, but for the cluster NGC 2210.

Fig. B.3 Same as in Fig. 2, but for the cluster NGC 2257.

Appendix C Photometric uncertainties as a function of radius

As shown in Figure C.1, due to crowding, the photometric errors increase sharply as the radius decreases within about 38 arcsec from the center of the cluster. Because of this, we select only stars within the range of 10 to 30 pc for the analysis in this work. We evaluate the photometric uncertainties near the MS by calculating the absolute value of the difference between the measured magnitudes of detected ASs and their input magnitudes. We calculate the corresponding color errors (as shown on the vertical axis of Figure A1) by adding the squared errors of the two passbands and then taking the square root.

Fig. C.1 Color uncertainties as a function of radius from the cluster center This figure depicts color uncertainties plotted against the radius from the center of the cluster. The shaded region highlights the radius range where color uncertainties remain approximately constant.

Appendix D Extinction residuals

Figure D.1 shows an example we used to evaluate the differential extinction. The top left panel shows an input differential extinction spatial distribution map based on the ASs (with colors representing their extinction values), while the top right panel shows the derived differential extinction map using the statistical methods described in Section 2.5. The difference between these two maps gives us the extinction residuals left by our method. Based on this approach, we evaluate the most likely color residuals when obtaining different differential extinctions, and these results are subsequently used to correct our ASs photometry.

Fig. D.1 Deriving differential reddening residuals. The upper panels show the distribution of simulated stars with differential extinctions (upper-left), and with their differential extinction distribution derived from statistical methods (upper-right, see Section 2.5). The lower panel displays the correlation between the color residuals (after we corrected the differential extinction) and the values of differential extinctions.

References

All Tables

All Figures

	Fig. 1 Stellar Number Density Profiles of Hodge 11, NGC 1841, NGC 2210, and NGC 2257. Best-fitting King profile shown by solid red curves. Solid blue line, dashed orange line, and dash-dotted green line indicate core radius, half-number radius, and tidal radius, respectively.
In the text

	Fig. 2 CMD of Hodge 11. Left: CMD shows all stars observed in the field of cluster Hodge 11. Middle: CMD displays stars classified as cluster members. Right: CMD features stars identified as field stars. In each panel, inset plots illustrate differences in coordinates between two observations taken approximately 7 years apart (top) and the spatial distribution of the corresponding stars.
In the text

Fig. 3 χ2 contour maps as a function of internal He spread, δY, and fraction of He-enriched stars, f2P . χ2 values are coded by colors. Best-fitting points (minimum χ2) are represented by red stars. Black solid lines indicate 1σ uncertainties, calculated from χ2 variations. The adopted model assumes a uniform He distribution. Panels a, b, c, and d correspond to clusters Hodge 11, NGC 1841, NGC 2210, and NGC 2257, respectively.

In the text

	Fig. 4 Same as in Fig. 3, with the adopted model assumes a bimodal He distribution.
In the text

	Fig. 5 Comparison of MS observations and synthetic MS. From left to right are comparisons between the observed MS, and the synthetic MS of the SSP, of the best-fitting MPs with a uniform He distribution, and of the best-fitting MPs with a bimodal He distribution. From top to bottom are clusters Hodge 11, NGC 1841, NGC 2210 and NGC 2257.
In the text

	Fig. 6 Comparison of observed and simulated HB. From left to right are the observed and simulated HBs for clusters Hodge 11, NGC 1841, NGC 2210 and NGC 2257. color bars indicate the helium mass fraction, Y, of the simulated HB stars.
In the text

Fig. 7 Internal helium spread and cluster mass relationships. Top: the internal helium spreads and the clusters present-day masses diagram, δY − log( M/M⊙). Red pentagrams are LMC clusters. Small grey dots are Milky Way GCs. Dark/light grey circles are SMC clusters. Bottom: the correlation between the internal helium spreads and the clusters initial masses, for Milky Way GCs, LMC and SMC clusters. The arrow indicates that we consider this value to represent only an upper limit. The arrow indicates that we consider this value to represent only an upper limit.

In the text

	Fig. B.1 Same as in Fig. 2, but for the cluster NGC 1841.
In the text

	Fig. B.2 Same as in Fig. 2, but for the cluster NGC 2210.
In the text

	Fig. B.3 Same as in Fig. 2, but for the cluster NGC 2257.
In the text

	Fig. C.1 Color uncertainties as a function of radius from the cluster center This figure depicts color uncertainties plotted against the radius from the center of the cluster. The shaded region highlights the radius range where color uncertainties remain approximately constant.
In the text

Fig. D.1 Deriving differential reddening residuals. The upper panels show the distribution of simulated stars with differential extinctions (upper-left), and with their differential extinction distribution derived from statistical methods (upper-right, see Section 2.5). The lower panel displays the correlation between the color residuals (after we corrected the differential extinction) and the values of differential extinctions.

In the text