© M. Tailo et al. 2022

1. Introduction

The search for solar-like oscillations in globular clusters (GCs) is difficult. Among the Galactic GC, M4 is the closest and the only one currently accessible to asteroseismic inference. Early attempts with ground-based instruments by Frandsen et al. (2007) reported hints of detection, but the low signal-to-noise ratio and the complexity of the spectra made the detection ambiguous. Similarly for another metal poor GC, NGC 6397, Stello et al. (2009) tried to detect solar-like oscillation but, also in this case, their results were uncertain due to low a signal-to-noise ratio.

The NASA K2 mission (Howell et al. 2014) observed M4 during its second observational campaign, providing ∼80 d of nearly continuous, high-precision photometric monitoring. Using K2 data, Miglio et al. (2016) reported a clear detection of solar-like oscillations in K giants in M4, obtaining estimates of mass, radius, and age compatible with estimates from other methods such as studies of eclipsing binaries (e.g. by Kaluzny et al. 2013), providing a crucial test of asteroseismology in the low-metallicity domain ([Fe/H] = −1.1, Marino et al. 2008; Carretta et al. 2013). However, given the limited sample of stars with detected oscillations (eight targets), Miglio et al. (2016) could not address key questions about M4, such as providing a robust estimate of the mass loss along the red giant branch (RGB) or investigating the features of the horizontal branch (HB) stars.

The crowding of a GC field and the drift of the spacecraft make the detection of solar-like oscillations challenging. Wallace et al. (2019) developed a procedure to mitigate these issues and extracted more than 4000 clean, high-quality time series of M4 stars, which are publicly available. In their analysis, Wallace et al. (2019) identified 55 stars – the majority on the RGB – showing potential evidence for solar-like oscillations. These high quality observations can help solve some of the questions left open by Miglio et al. (2016).

Being an old GC, M4 hosts multiple populations (Renzini et al. 2015; Bastian & Lardo 2018; Gratton et al. 2019) with different light elements abundance patterns and helium enhancement. In the case of M4, only two major groups are present: (i) a first generation of stars with light elements abundances compatible with field stars of a similar metallicity (e.g. high [O/Fe] and low [C+N/Fe] and [Na/Fe], with the standard helium abundance for its metallicity), and (ii) a second generation of stars with different chemical patterns, for example lower [O/Fe] and higher [C+N/Fe] and [Na/Fe] with a helium enhancement of 0.013 (Marino et al. 2008, 2011; Tailo et al. 2019).

In this Letter, we extend the work presented in Miglio et al. (2016) to the candidates identified in Wallace et al. (2019) with four objectives: (i) infer masses and radii of the stars in M4 and compare them with independent measurements; (ii) test the asteroseismic properties of metal poor stars; (iii) analyse the multiple populations in M4 (see Tailo et al. 2019) with asteroseismology; and (iv) estimate the integrated mass loss along the RGB by estimating the stellar parameters of red HB stars.

2. Data sets and observational constraints

We computed a power spectral density (PSD) for each star for which light curves from Wallace et al. (2019) were available. This procedure was carried out for a total of 4554 stars by means of the KADACS libraries (e.g. García et al. 2011). We performed a first skimming of the potential candidates by visually inspecting the PSDs and checking their location in the colour-magnitude diagram (CMD). Our final sample includes 54 stellar candidates located on the RGB and in the red part of the HB (rHB). We then fitted each PSD by means of the background modelling technique presented by Corsaro & De Ridder (2014), Corsaro et al. (2015, 2017) based on the DIAMONDS code. The background model adopted for each star was selected through a Bayesian model comparison process.

In addition to fitting the background, we obtain the position of the Gaussian envelope of the stellar oscillations and thus the frequency of maximum power, ν_max, for each star, while the large frequency separation, Δν, was estimated from the auto-correlation function (ACF) over the range of the PSD that contains the stellar oscillations. We have confirmed detection in 37 stars, six of which are located in the rHB. Three examples of the PSD fit are shown in the left part of Fig. 1. More details on the procedure are reported in Appendix A. The determined global asteroseismic parameters (ν_max and Δν) are listed in Table 1.

Fig. 1. Summary of the targets and spectra analysed in this work. Left: background fits (red curve) overlaid on the stellar PSD (grey) for the stars W837 (lower RGB), W1225 (HB), and W2665 (upper RGB). A smoothing of the PSD using a boxcar with width set to Δν is shown for reference as the black curve. The dot-dashed blue curves depict the Harvey-like components, while the horizontal dot-dashed yellow line and the dot-dashed magenta curve represent the level of the white noise and the Gaussian envelope of the oscillations, respectively. The dashed green line on top of the background fit shows the resulting fit when incorporating the oscillation envelope. We note here that the number of components in the model used for the background fit has been identified as part of the procedure (see Sect. 2 and Appendix A). Right: optical CMD of M4 from Stetson et al. (2019) data in the B and V bands. We show the position of the targets with good photometry and confirmed oscillation excess. The three stars in the left panels are identified as labelled.

We used stellar astrometry from Gaia and photometry in the U, B, V, and I bands from Stetson et al. (2019) to better characterise our targets. We identified cluster members by using the method by Cordoni et al. (2018), which is based on Gaia eDR3 parallaxes and proper motions (Gaia Collaboration 2021). Photometry has been corrected for differential reddening as in Milone et al. (2012). In the right panel of Fig. 1, we plotted the CMD in B − V versus V bands in which we highlight the selected targets, each of them being coloured according to its ν_max value. Our sample significantly extends, in number and ν_max range, the sample in Miglio et al. (2016). Five stars in our sample have uncertain photometry¹, and therefore we have removed them from our analysis.

We have exploited the capabilities of the C_U, B, I = (U − B)−(B − I) index (Milone et al. 2013; Monelli et al. 2013) to separate the two stellar populations in this cluster. In Table 1, we report the values of the corrected V magnitude and the value of the C_U, B, I index for all our targets with good photometry. The correct identification of the generation to which our RGB stars belong is necessary to correctly assess the integrated mass loss for the stars in this GC because our rHB stars all belong to the first generation (see Marino et al. 2011, and the first part of Appendix B). After the selection procedure, our final sample contained 26 RGB stars (16 s generation or 2G, and ten first generation or 1G) and six rHB stars.

A set of temperatures was derived from V − I colours as these are less affected by the shift in chemistry between the stellar populations in M4. In order to explore possible systematic uncertainties introduced by the choice of bands, we derived a second set of temperatures from the B − V colour. We adopted E(B − V) = 0.37 and E(V − I) = 0.53 from Hendricks et al. (2012). To obtain the values of T_eff for both combinations of bands, we used the bolometric corrections and the colour–T_eff relation from Casagrande & VandenBerg (2014), adopting [Fe/H] = −1.1 and [α/Fe] = +0.4 Marino et al. (2008). An iterative procedure was adopted to refine the results by using the asteroseismic surface gravity obtained from ν_max and T_eff. For consistency with the adopted colour–T_eff relations, we used the distance modulus (m − M)₀ = 11.20  ±  0.1 from Casagrande & VandenBerg (2014)² which – combined with the magnitudes of our stars – provide us with their luminosity and hence, by using the derived effective temperature, their photospheric radius (R_CMD).

3. Masses and radii

Combining the above decided asteroseismic and non-asteroseismic parameters, we estimated stellar masses and radii with classical scaling relations. Combining Δν ∝ ρ^1/2 and (see e.g. Brown et al. 1991; Frandsen et al. 2007; Chaplin & Miglio 2013; Miglio et al. 2012, 2016) with , we obtain the set of equations reported below:

(1)

(2)

(3)

(4)

We would like to note that Δν is not taken at nominal value: multiplicative corrections to the relation between Δν and ρ^1/2 (f_Δν) were applied following the approach described in Rodrigues et al. (2017) which, so far, has yielded masses and radii showing no systematic deviations to within a few percent of independent estimates (see e.g. Miglio et al. 2016; Rodrigues et al. 2017; Handberg et al. 2017; Brogaard et al. 2018, who partially revisited the work by Gaulme et al. 2016; Buldgen et al. 2019). Our computations of f_Δν are described in Appendix C, while the full set of parameters is reported in Table D.1.

We may also estimate the stellar radii with the following:

(5)

For solar reference values, we adopted Δν_⊙ = 135.1 μHz, ν_max, ⊙ = 3090 μHz, and T_eff, ⊙ = 5777 K (Huber et al. 2013). Finally the 1σ interval on each quantity has been obtained with standard error propagation methods, starting from the 1σ uncertainty of the parameters involved in each equation.

As an independent check for potential biases in the seismically inferred masses and radii, we compared the asteroseismic radius (R_S) and the photometric one of all of our targets. For RGB stars, we find that ⟨(R_S − R_CMD)/σ⟩ = 0.11 ± 0.15 of the combined error. This lends confidence to both the mass and distance scale, given the strong correlation between seismic masses and radii (see e.g. Basu et al. 2011; Wu et al. 2014; Khan et al. 2019; Miglio et al. 2021, and references therein). For rHB stars, however, the situation is different as we find significant differences in radii. For these stars we have ⟨(R_S − R_CMD)/σ⟩ = 1.58 ± 0.77 and at the same time the mass of our HB stars calculated with the scalings involving Δν (Eqs. (1), (2), and (4)) do not agree with each other. A thorough exploration of possible systematics (see Appendix F) points to an issue related to measuring Δν in rHB.

We checked the robustness of our Δν measurement via a simple test with the Kepler data. We degraded the entire 4 yr duration of the Kepler observations for two rHB stars – namely KIC 8694070 and KIC 11450315 – by diving the full time series into 80-day chunks and we derived Δν from these shorter data sets. We find that the measured Δν fluctuates up to 40%, leading one to overestimate (or underestimate) the mass and the radius of the HB stars by up to ∼8 and ∼3 times, respectively. We give more details on this in Appendix A.1. We believe that this is due to the intrinsic complexity of the HB oscillation spectra (Matteuzzi et al., in prep.), which combined with the short duration of K2 observations make accurate measurements of Δν difficult. This kind of instability is not present in RGB stars which, on the other hand, present no issues when measuring Δν. Therefore, when discussing rHB stars, we use Eq. (3), which does not include Δν and this unstable behaviour, to estimate the stellar parameters of rHB stars. A thorough investigation of systematic effects (in T_eff, distance, Δν, and ν_max) is presented in Appendix F.

3.1. Mass of the RGB stars in M4

We report the mass of our RGB targets in the four panels of Fig. 2 as the collection of red and blue points for 1G and 2G stars, respectively (see Appendix B for more details). The weighted average mass for our 26 RGB stars (⟨M₁⟩ = 0.842 ± 0.035 M_⊙, ⟨M₂⟩ = 0.847 ± 0.015 M_⊙, ⟨M₃⟩ = 0.867 ± 0.007 M_⊙, and ⟨M₄⟩ = 0.863 ± 0.024 M_⊙) are shown as the black dashed lines in each panel. Combining these values with the isochrones from Tailo (2016) and Tailo et al. (2020), we found an age in the range 11−12 Gyr, with the exact value depending on the adopted equation for the RGB mass³. These four values of mass and age agree with other published values for this cluster (Dotter et al. 2010; Kaluzny et al. 2013; VandenBerg et al. 2013; Miglio et al. 2016; Tailo et al. 2019; Jang et al. 2019). It is important to note that the different expressions linking the seismic parameters with the classical ones involve the value of Δν, ν_max, T_eff, and luminosity (L) with different powers. This is the main reason for the difference in the precision with which the stellar parameters were derived with the different scaling relations. Within the current uncertainties, no mass variation is observed over the magnitude range examined.

Fig. 2. Mass of our targets estimated with Eqs. (1)–(4), in descending order. The black dashed lines represent the average mass of the RGB population, while the coloured ones refer to the same quantity for the corresponding groups.

We also estimated the mean mass of 1G and 2G RGB stars; however, we do not find a significant mass difference given the current systematic and statistical uncertainties (see Appendices B and F for more details). This is not surprising as the predicted theoretical mass difference between the two populations is estimated to be 0.017 M_⊙ (Tailo et al. 2019), with the 1G being composed of more massive stars. This mass difference has been obtained assuming that the two populations are strictly coeval; therefore, it originates only from the effect of chemical enrichment on the evolution timescale. However, all of the most successful scenarios for the formation of multiple stellar populations predict that the two generations of stars form within a small interval of time, compared to the age of the cluster. Such an interval should be a few 10⁸ yr, at the largest, depending on the scenario (see Renzini et al. 2015; D’Antona et al. 2016; Bastian & Lardo 2018; Kim & Lee 2018). Therefore, when one takes into consideration that the 2G is the youngest, the mass difference reduces. Our results seem to suggest this; however, our uncertainties prevent us from saying more.

3.2. The mass of the HB stars

The masses for our six HB stars are reported as the orange diamonds in the third panel of Fig. 2. The average value is ⟨M_HB⟩ = 0.632 ± 0.024 M_⊙, which is in good agreement with other estimates of the average rHB mass for this cluster (Tailo et al. 2019; Jang et al. 2019). Even if the number of targets is small, the scatter shown by these values does not seem compatible with a single mass. This is shown in the lower panel of Fig. 3, where we compare the Gaussian kernel-density distribution of our stars (orange) with a single mass distribution (dashed black), centred at the ⟨M_HB⟩ with the corresponding σ.

Fig. 3. Properties of the HB stars in our sample. Top: temperature of our rHB targets as a function of their mass. Bottom: Gaussian kernel distribution for the data (orange). The black dashed line represents a single-mass population centred at the value of ⟨MHB⟩ with its corresponding sigma.

An inspection of the properties of these stars reveals the presence of a clear trend between mass and effective temperature (top panel of Fig. 3), with an estimated Spearman r index, r ∼ −0.94 (p < 0.005), indicating that this gradient is significant at the ∼3σ level. We note that this behaviour is consistent with theoretical models, which predict hotter stars for smaller masses (with a smaller H envelope mass, see Girardi 2016, for instance), and, in such a case, this would be the first direct observation of HB extension in terms of mass variations.

3.3. Integrated mass loss

The average integrated RGB mass loss can be obtained by subtracting the average mass of the rHB stars from that of the RGB ones. To get the correct value of integrated mass loss, the rHB stars have to be compared with their progenitor in the 1G, along the RGB (with an average mass of 0.859 M_⊙, see Appendix B). The estimated average integrated mass loss of the 1G stars is therefore ΔM = 0.227 ± 0.028 M_⊙, which is compatible with the one in Tailo et al. (2019; ΔM = 0.209  ±  0.024 M_⊙). If this value is converted into a mass loss rate (as in Fig. 16 in Tailo et al. 2020), the associated Reimers (1975) mass loss parameter is η_R ∼ 0.48, while using Padova stellar models (Bressan et al. 2012) a value closer to 0.50 is obtained.

Adding this result to the previous determinations of the mass loss rate in stellar clusters using asteroseismology (Miglio et al. 2012, 2021; Stello et al. 2016; Handberg et al. 2017) allows us to cover a wide range in age and metallicity. In Fig. 4, we show the mass loss expected from stellar models from the Padova database (Bressan et al. 2012), with a different chemical composition, age, and two values of the Reimers (1975) mass loss parameter η_R = 0.25; 0.50. These predictions were compared with the integrated mass loss estimated from asteroseismic studies of stellar clusters (including the value for M4 derived in this Letter) and of the high-[α/Fe] red giants in the Kepler field (Miglio et al. 2021). The agreement between the asteroseismic estimates of mass loss and the estimate coming from fitting the photometric HB data with models of stellar populations (Tailo et al. 2019; Jang et al. 2019) corroborates the hypothesis from Tailo et al. (2020, 2021) that mass loss in old GCs is a substantially different phenomenon than in open clusters and in field stars which, if one assumes Reimers’ parameterisation, is properly described only invoking a higher mass loss parameter.

Fig. 4. Mass loss measured using asteroseismic observables for M4, NGC 6791 (Miglio et al. 2012), NGC 6819 (Miglio et al. 2012; Handberg et al. 2017, H17), M 67 (Stello et al. 2016), and the α-enhanced stars in the Kepler field (Kα, Miglio et al. 2021) as function of age (left) and metallicity (right). Two sets of models for different (Reimers 1975) mass loss parameter, ηR = 0.25, 0.50, are included for comparison purposes.

This divergent behaviour could be connected to different environmental factors (an interaction with the Galaxy or a different formation environment), be the footprint of an early dynamical interaction between the stars in GCs (stellar collisions and binaries), be the product of some still poorly understood physics inside RGB stars (such as the mechanism proposed in Fusi-Pecci & Renzini 1975, 1976), or even a combination of all the previous. Indeed, the extensive study performed in Tailo et al. (2020, 2021) shows that the large population of GCs analysed lies on a relation systematically higher than the one for the open clusters and the field stars of a comparable metallicity (where the comparison is possible). Furthermore, said relation is compatible with the ones obtained for dwarf spheroidal galaxies (e.g. Salaris et al. 2013; Savino et al. 2019). This suggest that RGB mass loss in old stellar association is somewhat universal.

4. Conclusions

This study presents a new asteroseismic analysis of targets located in M4 which features the largest sample in these kinds of studies for GCs (37 stars). We combined asteroseismic measurements of ν_max and Δν with the Stetson et al. (2019) photometric catalogue to provide a direct measurement of 32 stars, six of which are rHB stars.

Our results confirm the validity of asteroseismology in the low-metallicity regime, thus giving more solidity to the results in Valentini et al. (2019), Montalbán et al. (2021), and Miglio et al. (2016, 2021), and providing a direct estimate of the RGB mass loss in this cluster which is in agreement with the estimate derived from the fitting of the HB photometric data with stellar populations models (Tailo et al. 2019; Jang et al. 2019). Although we separated the stellar populations in this cluster efficiently, our uncertainties prevent us from commenting on the mass difference between the two stellar populations. However, this work paves the way for future studies using a more detailed asteroseismic analysis of the individual mode frequencies which will allow for more precise mass estimates (see e.g. Montalbán et al. 2021) that are less reliant on the method used to obtain T_eff. Finally, the gradient we observe in our HB stars is the first, direct, and model-independent observation showing the colour-temperature-mass correlation known to exists in HB stars from the theoretical models.

Acknowledgments

M.T., A.M., J.M., A.S., G.C. acknowledge support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, https://www.asterochronometry.eu, G.A. No. 772293). A.P.M. acknowledges support from the European Union’s Horizon 2020 research innovation programme (Grant Agreement ERC-StG 2016, No. 716082 ‘GALFOR’, PI: Milone). Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant agreement No. DNRF106). E.C. acknowledges support from PLATO ASI-INAF agreement No. 2015-019-R.1-2018.

References

Appendix A: Obtaining the global asteroseismic parameters

To assess the presence of an oscillation power excess, we computed a PSD for each star for which light curves from Wallace et al. (2019) were available. This procedure was carried out for a total of 4554 stars by means of the KADACS libraries (e.g. García et al. 2011). A skimming of potential candidates was done by a visual inspection of the PSDs in conjunction with the location of the star in the CMD. We ended up with a list of 54 stellar candidates, with evolutionary stages spanning from the low part of the RGB to the HB.

We proceeded by performing an actual fit to each PSD of the candidate stars by means of the background modelling technique presented by Corsaro & De Ridder (2014), Corsaro et al. (2015), and Corsaro et al. (2017), which involves the use of the public code DIAMONDS+BACKGROUND based on the Bayesian nested sampling Monte Carlo algorithm (Skilling 2004). The background model adopted for each star was selected through a Bayesian model comparison process that uses Bayesian evidence as a statistical quantity for model selection. Bayesian evidence is a direct output of a statistical inference performed by means of the DIAMONDS code. As part of its output, the Bayesian procedure also indicates the number of Harvey-like components (Harvey et al. 1985) needed to fit the background of the PSD. We have found that either a model consisting of a single or a double Harvey-like component represents an adequate choice for the data sets provided by Wallace et al. (2019). In this case, these Harvey-like components mostly refer to the presence of a granulation-activity signal. In addition to the Harvey-like components, we have included a Gaussian function to mimic the presence of a power excess caused by stellar oscillations. The location of the power excess was initialised by means of an input guess for the frequency of maximum oscillation power, ν_max, which was extrapolated by combining information from a visual inspection of the PSD and the position of the star in the CMD, hence its log g. An example of the fit produced is shown in Fig. 1 for stars in three different locations of the CMD. Here we can observe how the number of significant Harvey-like components decreases from two to one as the oscillation power excess moves to lower frequencies. The limited frequency resolution of the datasets increasingly penalises our capability to infer the background properties because the usable frequency range to fit a single Harvey-like profile shrinks as the star evolves along the RGB.

The presence of an oscillation power excess for the stars that were considered as potential candidates was further confirmed through the adoption of a Bayesian model comparison, which was once again performed by means of the DIAMONDS code. This procedure was carried out by comparing the Bayesian evidence of a model that incorporates the Gaussian envelope in addition to the Harvey-like functions, and of another model that only contains the latter ones. We considered the minimum for a detection in a single star to be a weak evidence condition (Trotta 2008), corresponding to lnB₀₁ ≃ 1.0, where lnB₀₁ is the natural logarithm of the Bayes’ factor, that is to say the ratio between the Bayesian evidence of the models with and without the oscillation envelope. We have detected solar-like oscillations in 37 stars, with evolutionary stages spanning from the RGB to the HB.

A.1. Measuring Δν

In addition to testing the detection for solar-like oscillations, the fitting of the background allowed us to obtain the value of ν_max for each star. However, the estimation of the large frequency separation Δν required us to compute an ACF as well over the range of the PSD that contains the stellar oscillations. This procedure was carried out by adopting a search range in Δν that is centred around an input guess based on ν_max through the Δν–ν_max relation by Stello et al. (2009), and having an extent of about ±50 % with respect to the input guess. In most cases, the ACF produces a clear peak that is relatively close to the input guess from the empirical relation. For completeness, in Fig. A.1 we compare the ν_max–Δν relation for our targets with the one for the Kepler stars. In the plot, the latter stars are colour-coded according to their mass from Yu et al. (2018). We also included the stars from the open clusters M67 (Stello et al. 2016) and NGC 6819 (Handberg et al. 2017).

Fig. A.1. νmax–Δν relation for our targets divided by groups. We compare our targets with the Kepler stars, which are colour-coded according their mass (Yu et al. 2018), and the stars in NGC 6819 from Handberg et al. (2017, H17) and M67 from Stello et al. (2016).

However, for some of the stars belonging to the HB, the obtained Δν value may deviate significantly (even up to about 40 %) from what one would expect. Since HB stars have more complex oscillation spectra than their RGB counterparts, which causes a less clear Δν signal in the ACF, we decided to evaluate whether these deviations can be affected by systematics related to the relatively short time span of the observations. For this purpose, we conducted a test involving two red HB stars that were observed by NASA Kepler for more than 4 years and for which a reliable estimate of Δν can be obtained, namely KIC 8694070, and KIC 11450315. We created multiple realisations of their PSDs by dividing their light curves into chunks with 80 d length each to simulate the observations from Wallace et al. (2019). We obtained a total of 18 chunks for KIC 8694070 and 16 chunks for KIC 11450315, as shown in Fig. A.2. We then fitted the background signal in the PSD of each chunk following the same approach presented in Appendix A, hence we measured Δν accordingly. The results plotted in Fig. A.2 show that the individual Δν measurements suffer from large scatter compared to the reference value, which was obtained from the full-length light curve; although, on average, this behaviour is not significant when the sample dispersion is taken into account. Most important, however, is that the underestimation (overestimation) can be up to about 40 % (20 %) of the reference Δν. This in turn implies that the corresponding asteroseismic radius may be overestimated (underestimated) by up to a factor of about three, which is in qualitative and quantitative agreement with what we obtain by comparison to the CMD estimates of the stellar radius for some of our cluster HB stars.

Fig. A.2. Deviations in the percentage of Δν for individual chunks each corresponding to 80 day-long time series of the core-He-burning stars KIC 8694070 (top) and KIC 11450315 (bottom) observed by Kepler. The reference value of Δν, Δν0, is obtained from the full-length light curve. Measurements of Δν obtained from individual chunks are shown by orange dots, with associated 1-σ error bars. The horizontal dashed line corresponds to Δν0, with its 1-σ confidence region marked by the teal shading. The grey shading indicates the 1-σ confidence region of the sample of chunk measurements.

In addition to this issue relating to the HB stars, we find the possibility of a second effect that intervenes when measuring Δν of the RGB stars worth mentioning. In fact there could be a small systematic shift in the value obtained with ACF compared to the one obtained by fitting individual frequencies (when the accuracy is high enough Viani et al. 2019). According to Fig. 1 in Khan et al. (2019), this shift is ≃ + 1%, therefore falling into the possible systematic effects discussed in Fig. F.1.

Appendix B: Multiple populations of stars in M4

The main feature of GCs, such as M4, is that they host a collection of different stellar populations, differing in light elements and helium abundances (Renzini et al. 2015; Bastian & Lardo 2018; Gratton et al. 2019). Generally the populations in a GC can be divided into two groups: a first generation (1G), having light element abundances compatible with the ones of the field stars of a similar metallicity (high [O/Fe] and [Mg/Fe], low [C+N/Fe], [Na/Fe], [Si/Fe], and [Al/Fe], with the standard helium abundance for its metallicity), and a second generation (2G, generally divided into other sub-groups, see Carretta et al. 2009; Milone et al. 2015, 2017; Bastian & Lardo 2018; Milone et al. 2018; Marino et al. 2019; Gratton et al. 2019, and references therein)showing radically different chemical patterns (lower [O/Fe] and [Mg/Fe], higher [C+N/Fe], [Na/Fe], and [Al/Fe], with helium enhancement that can reach 0.35 in some cases, see an example for NGC 2808 in Fig. 2 of Gratton et al. 2019, and references therein). It is therefore natural to check what populations our stars belong to, and if we are able to detect a significant mass difference between the populations.

Among our targets, the most straightforward to identify are the six HB stars. Being located in the red part of the HB (see Fig. 1), they can be safely identified as 1G stars following the results of Marino et al. (2011, who found all rHB in this cluster to have high [O/Fe] and low [Na/Fe]). This can be further reinforced by cross-matching our sample with the one in Marino et al. (2011): we find two stars in common (W1068 and W2887) which have [Na/Fe] = −0.09 and −0.06, respectively. A value of [Na/Fe] close to zero indicates that the star belongs to the 1G (see e.g. Gratton et al. 2001; Marino et al. 2008; Carretta et al. 2009), as 2G stars tend to have enhanced abundances of [Na/Fe] (and corresponding depleted abundances of [O/Fe]). This corroborates the identification of all rHB stars in our sample as 1G stars.

To separate the RGB stars between 1G and 2G, we exploited the capability of the C_U, B, I = (U − B)−(B − I) index (Milone et al. 2013; Monelli et al. 2013). We summarise the procedure in Figure B.1. To summarise, we applied a Gaussian mixture model (GMM) to the verticalised RGB in the C_U, B, I versus V diagram (left and upper right panels of Figure B.1). A free parameter for the GMM is the number of clusters it finds. We tested the number of components across a range of two to 15 and based on an evaluation of the different Bayesian information criteria values that three clusters made up the model that best described the data in this dimension. Our best-fit model has three components: two representing the main stellar populations in the cluster 1G (red) and 2G (blue) and one⁴ (green) identifying the few binaries in our catalogue. This is in agreement with most literature sources that identify two stellar populations in this GC (see Marino et al. 2008; Carretta et al. 2009, 2013; Marino et al. 2011; Milone et al. 2017; Marino et al. 2017; Lardo et al. 2017; Milone et al. 2018; Tailo et al. 2019, and references therein).

Fig. B.1. Summary of the method used to separate the two stellar populations along the RGB of M4. Left panel: The photometric data in the CU, B, I vs. V plane. Our targets are highlighted as orange points. The two dashed lines mark the location of the 10th and 90th percentiles of the RGB stars for the left and the right line, respectively. Upper right panel: The verticalised RGB (ΔCU, B, I vs. V) of M4. The points are divided between 1G (red) and 2G (blue) according to the GMM reported in the bottom right panel. The green points, seen far from the main RGB, represent the binary population. Bottom right panel: The GMM model used to identify the two stellar populations. The two most prominent Gaussian curves are the ones for the 1G and 2G stars.

Our targets can be divided as follows: ten 1G stars and 16 2G stars. We took special care when verifying that the identification of our targets was not significantly affected by random factors. In order to assess whether the identification was wildly affected by the randomisation inert in the GMM computation, we tested whether the identification was robust across 20 different initialisations, finding that none of our targets switched identification. We highlight the position of our targets in the panels as large dots, colour-coded per the population.

We obtained a further confirmation of the separation between the 1G and the 2G stars in our target sample by cross-matching our stars with the ones in Marino et al. (2008) and Carretta et al. (2009). Only 12 stars are in common between these spectroscopic works and our sample and the spectroscopic and photometric identification agree very well, as reported in Table 1. One star (W2022) is classified as 1G in Marino et al. (2008) and 2G in Carretta et al. (2009); we agree with the former, which is based on higher resolution spectra. We note that the ratio of 1G to 2G stars (∼40% and ∼60%, respectively) is compatible to the ratio found in the literature by Marino et al. (2008), Carretta et al. (2009), Carretta et al. (2013), Marino et al. (2011), Milone et al. (2017), Marino et al. (2017), and Lardo et al. (2017).

B.1. The average mass of the multiple populations

With our targets separated into first and second generation stars, we re-evaluated the average mass values. We found the average mass values reported in Table B.1. We report the two separate sets of average mass values in the panels in Fig. 2 as the red and blue solid lines for the 1G and the 2G, respectively, together with their 1σ interval. Considering our uncertainties on the derived masses, which are both statistical and systematical, our results do not allow one to identify significant mass differences between the two populations (i.e. the mass difference is comparable with the errors). This is, however, unsurprising. Due to the low helium enhancement of the two stellar populations in M4 (0.013) and the low extension of the [Na/Fe]-[O/Fe], [(C+N)/Fe]-[O/Fe], and [Mg/Fe]-[Al/Fe] anti-correlations’ diagrams, compared to more massive and complex clusters such as NGC 2808 or M13, the difference between the evolving mass of the 1G and the 2G is predicted to be ≤0.017 (Tailo et al. 2019); this sensitivity is hard to reach with few targets in both populations when the masses are derived purely using the asteroseismic scaling relations instead of the more constraining individual mode frequencies.

Appendix C: Evaluation of correction factor f_Δν

Corrections to scaling relations linking Δν to the mean density have been obtained with the iterative procedure described in the following. For the RGB stars, we calculated two grids of stellar models with the MESA software (Paxton et al. 2011, 2013, 2015, 2018, 2019). Each grid was built with a Cartesian approach and spans the following dimensions as follows: M from 0.50 to 1.20 M_⊙ with a step of 0.02M_⊙; [Fe/H] = −1.1, [α/Fe] = +0.4, and standard helium (Y∼0.25) for the 1G and M from 0.50 to 1.20 M_⊙ with a step of 0.02M_⊙; and [Fe/H] = −1.1, [α/Fe] = +0.3, and Y = 0.263 for the 2G. The increased helium abundance and the lower [α/Fe] were needed to simulate the different chemical abundance patterns observed in the 2G stars (Marino et al. 2008; Carretta et al. 2009, 2013; Tailo et al. 2019). Each model was evolved up to the helium flash.

The grid for the HB stars covers a range of initial masses that ranges from M = 0.82 M_⊙ to 0.90 M_⊙ with a step of 0.01M_⊙ and the same [Fe/H] and [α/Fe] from the 1G grid. We assumed a mass loss rate corresponding to a Reimers (1975) parameter of η_R = 0.47, which is in line with the predictions of both Tailo et al. (2020) and McDonald & Zijlstra (2015). The HB models were evolved up to the end of the HB phase, including the helium flash, and they arrived at the start of the HB phase with a mass in the M = 0.56 M_⊙ to 0.72 M_⊙ range.

We then computed the individual oscillation frequencies for each model in both grids using the GYRE software (Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020) to calculate individual oscillation frequencies. The correction to the Δν for each track was then obtained with the method from Miglio et al. (2016) and Rodrigues et al. (2017). Once each single track had been built this way, we used the mass calculated from Eq. 3 and the ν_max of each star to locate it along the track. For the RGB stars, the results of this iterative procedure are represented in the upper panel of Fig. C.1. To calculate the correction to be applied to the rHB stars, we located our stars in the track corresponding to the HB mass obtained from Eq. 3, and then located the position

Fig. C.1. Correction values for our RGB targets (left), divided by populations, and our rHB stars (right). The solid lines are the corresponding models from which the corrections are drawn.

of the tracks using the value of ν_max as well. The results for the rHB stars is plotted in the lower panel of Fig. C.1.

Compared with Miglio et al. (2016), we used a larger range of correction values. Miglio et al. (2016) and collaborators used a fixed value of f_Δν = 0.956 which is fine for the limited range of log(ν_max) of their sample. In this work, changes to their procedure are necessary given the greater range of ν_max and due to the luminosities we are exploring. On the other hand, the values we adopted for the rHB stars are quite different from the f_Δν, HB = 1.029 adopted by Miglio et al. (2016) and collaborators. This is due to the lower mass loss they employed in their models. However, we think that it is worth noting that models with lower mass loss – resulting in more massive HB stars – do not reach values of ν_max compatible with the ones of some of our stars. The entire range of f_Δν adopted is reported in Table D.1.

Appendix D: Full masses and radii sets

We report in Table D.1 the complete set of obtained masses, together with the photometric and asteroseismic radii for our sample of stars. Two sets of values are reported with uncorrected and corrected Δν, respectively.

Appendix E: Comparison of the radii

As shown in Appendix A.1 and in Fig. A.2, the value of Δν can vary wildly, especially for the HB stars. This causes the masses and radii of the stars in our sample to vary as well, leading to the severe underestimation or overestimation of both.

To conduct a simple yet effective quality check of our results, we compared the radius obtained from the asteroseismic parameters (R_S) with the one obtained from the photometry (R_CMD). The results are plotted in Fig. E.1 for the RGB (top) and the HB stars (bottom), respectively. We see that after the correction on Δν was applied, the radius of most of our RGB targets fell within |⟨(R_S − R_CMD)/σ⟩| < 1. The only exception is one 2G star (W4488) which is located at ⟨(R_S − R_CMD)/σ⟩∼ − 3σ. Since the mass obtained from Eq. 3 is 0.836 ± 0.072 M_⊙, which is in line with the average of its 2G counterparts, we have strong hints that the culprit of this large discrepancy is a wrong measurement of Δν. If we remove it from the sample, then the average mass of our RGB stars becomes 0.879 ± 0.025M_⊙, 0.836 ± 0.012M_⊙, 0.862 ± 0.008M_⊙, and 0.882 ± 0.026M_⊙ for Eq. 1 to 4, respectively.

Fig. E.1. Difference between the photometric (RCMD) and the asteroseismic (RS) radii for our RGB (top) and rHB (bottom) stars. The black dashed line in the upper panel is the average dispersion of the entire sample, while the coloured lines represent the average value of each population with their 1σ intervals as the shaded areas.

The results of this check for the rHB stars in our sample yields a different result, however. Because of the difficulties associated with measuring Δν, the radii of the rHB stars vary up to the |⟨(R_S − R_CMD)/σ⟩| < 5 level, with only two (namely W1225 and W2386) falling in the 1σ range. The points in the bottom panel of Fig. E.1 are also colour-coded according the value of M₁, showing that when the radius is overestimated the mass follows as well (and vice versa).

We carried out an additional test by searching for the value of f_Δν that gives us |⟨(R_S − R_CMD)/σ⟩| < 1 and average masses from Eq. 1 to 4 in agreement with each other for our six rHB stars. We found that with f_Δν ∼ 0.940, these two conditions are reasonably satisfied. Indeed, we obtained ⟨(R_S − R_CMD)/σ⟩= − 0.19 ± 0.76 and a set of average masses of 0.569 ± 0.298M_⊙, 0.512 ± 0.096M_⊙, 0.632 ± 0.024M_⊙, and 0.606 ± 0.168M_⊙. However, the correction needed is not compatible with the one predicted from the model plotted in Fig. C.1 and in other literature studies (Miglio et al. 2013; Rodrigues et al. 2017) for low-mass helium-burning stars. Finally, if we impose R_CMD = R_S, we obtain values for f_Δν that vary wildly and that are very different from the ones that can be obtained from the models in Fig. C.1. We can therefore exclude any major role of the correction to Δν in the large spread in ⟨(R_S − R_CMD)/σ⟩ observed among our rHB stars.

Based on these considerations, we only used the results following Eq. 3 from HB stars as they are not dependent on Δν. Regardless, the mass estimates obtained with all four equations are reported in Table D.1 for completeness.

Appendix F: Systematics

We discuss several systematical effects that can affect our results in this subsection. The results of this exploration are collectively reported in Figure F.1. Each row starting from the second one explores a different systematical shift and its effect on the average mass determination (first four columns) and on ⟨(R_S − R_CMD)/σ⟩ (last column) of the individual groups of stars in our sample⁵. The first row reports the non-shifted mass and scatter values that are also reported in each panel as the transparent points for clarity and ease of comparison. In each panel describing a mass shift, we report independent estimates coming from different sources as the coloured dashed line: the extrapolated mass from the eclipsing binaries (EBs) in Kaluzny et al. (2013, 0.85M_⊙, first column) and the mass of the 1G (red), the 2G (blue), and the rHB stars (orange) from Tailo et al. (2019, values are 0.85M_⊙, 0.833M_⊙, 0.64M_⊙ for the second to fourth columns, respectively) with their respective 1σ interval.

Fig. F.1. Systematic shift of the average mass and radii dispersion for different groups of stars. The first row reports the non-shifted values, which are also reported as transparent points in each panel. Each row past the first one explores the effects of a shift in a given parameter of its 1σ interval, as labelled. The last two rows report the non-corrected mass and dispersion average values, and what happens when the Temperatures coming from the B − V bands are used. The dashed lines represent the mass derived from the EB results in Kaluzny et al. (2013, 0.85M⊙, first column); the mass of the 1G, the 2G, and the HB from; Tailo et al. (2019, the values are 0.85M⊙, 0.833M⊙, and 0.64M⊙ for the second to fourth columns, respectively); and the shaded areas mark their 1σ intervals.

We explored what happens to our mass and ⟨(R_S − R_CMD)/σ⟩ estimates when the four main parameters in the scaling relations are shifted by ±1σ. We report the results in the second to fourth rows of Fig. F.1 for a shift in T_eff, distance modulus (hence luminosity), ν_max, and Δν, respectively. The variation of the average mass is proportional to the power that appears in the specific parameter for the scaling equations. This means that the mass shift in Δν, for example, is more pronounced for Eq. 1 than Eq. 2, and so on. Similar considerations can be drawn for ⟨(R_S − R_CMD)/σ⟩, where the shift for the RGB radii is significant, but does not go past the |⟨(R_S − R_CMD)/σ⟩| > 1 limit.

The values of the correction to Δν can be somewhat model-specific for several reasons connected to the exact physics used to calculate the underlying stellar models. To explore the systematic introduced by the procedure to obtain the correction, we shifted the values of f_Δν reported in Table D.1 by ±0.01. The results are reported in the third-to-last row of Fig. F.1. We see that, while a shift is present in all groups, in the majority of cases it is within 1σ. The same can be said for the shift in ⟨(R_S − R_CMD)/σ⟩ in the last column of the same row. For completeness, in the second-to-last row of the figure, we explore what happens when corrections are not taken into account (i.e. f_Δν = 1.0 for all stars).

When studying an old GC, which is made up of multiple stellar populations, differing in light elements abundance patterns, it is important to check whether the bands used in the temperature evaluation are affected by the star-to-star variation of the chemical patterns. To show this, we compared the temperatures estimates derived from the B − V and V − I colours. The results of this comparison are plotted in Fig. F.2. We see that the star-to-star variation of the temperature is well within the 1σ range; however, the shift in average temperature of the two stellar populations, while small – ∼27 K for the 1G and ∼55 K for the 2G – is significant and collectively affects the 2G stars more. For the 2G stars, this shift occurs as the B band is more sensitive to the changes in [C/Fe] and [N/Fe]. This small temperature difference is enough to alter the average mass difference of the two populations.

Fig. F.2. Shift in temperatures for our RGB targets when the estimate is carried out with the B − V or the V − I colours. The dashed lines are the average values for the two populations with their respective 1σ interval.

We show what happens to our average mass estimates when we adopted the temperatures coming from the B − V colours in the bottom row of Fig. F.1. We see that the agreement with the independent mass estimate for M4 gets slightly worse; this is because the lower temperature values produce slightly higher mass estimates. The shift is greater for the 2G RGB stars while within the 1σ range for the 1G RGB and rHB stars. Since the majority of stars are 2G stars in this case, they also affect the agreement of the combined sample. Similar considerations hold for the radii dispersion, shown in the rightmost panel of the bottom row in the Figure.

In the case of M4, studied in this work, the variations are minimal; however, in more complex and interesting cases, such as 47 Tuc (see e.g. Carretta et al. 2013; Milone et al. 2017, for a showcase), these effects can be relevant for the results. We then conclude that the best results are achieved with a combination of bands less affected by the abundance pattern shifts in the multiple populations.

All Tables

All Figures

Fig. 1. Summary of the targets and spectra analysed in this work. Left: background fits (red curve) overlaid on the stellar PSD (grey) for the stars W837 (lower RGB), W1225 (HB), and W2665 (upper RGB). A smoothing of the PSD using a boxcar with width set to Δν is shown for reference as the black curve. The dot-dashed blue curves depict the Harvey-like components, while the horizontal dot-dashed yellow line and the dot-dashed magenta curve represent the level of the white noise and the Gaussian envelope of the oscillations, respectively. The dashed green line on top of the background fit shows the resulting fit when incorporating the oscillation envelope. We note here that the number of components in the model used for the background fit has been identified as part of the procedure (see Sect. 2 and Appendix A). Right: optical CMD of M4 from Stetson et al. (2019) data in the B and V bands. We show the position of the targets with good photometry and confirmed oscillation excess. The three stars in the left panels are identified as labelled.

In the text

	Fig. 2. Mass of our targets estimated with Eqs. (1)–(4), in descending order. The black dashed lines represent the average mass of the RGB population, while the coloured ones refer to the same quantity for the corresponding groups.
In the text

	Fig. 3. Properties of the HB stars in our sample. Top: temperature of our rHB targets as a function of their mass. Bottom: Gaussian kernel distribution for the data (orange). The black dashed line represents a single-mass population centred at the value of ⟨MHB⟩ with its corresponding sigma.
In the text

Fig. 4. Mass loss measured using asteroseismic observables for M4, NGC 6791 (Miglio et al. 2012), NGC 6819 (Miglio et al. 2012; Handberg et al. 2017, H17), M 67 (Stello et al. 2016), and the α-enhanced stars in the Kepler field (Kα, Miglio et al. 2021) as function of age (left) and metallicity (right). Two sets of models for different (Reimers 1975) mass loss parameter, ηR = 0.25, 0.50, are included for comparison purposes.

In the text

	Fig. A.1. νmax–Δν relation for our targets divided by groups. We compare our targets with the Kepler stars, which are colour-coded according their mass (Yu et al. 2018), and the stars in NGC 6819 from Handberg et al. (2017, H17) and M67 from Stello et al. (2016).
In the text

Fig. A.2. Deviations in the percentage of Δν for individual chunks each corresponding to 80 day-long time series of the core-He-burning stars KIC 8694070 (top) and KIC 11450315 (bottom) observed by Kepler. The reference value of Δν, Δν0, is obtained from the full-length light curve. Measurements of Δν obtained from individual chunks are shown by orange dots, with associated 1-σ error bars. The horizontal dashed line corresponds to Δν0, with its 1-σ confidence region marked by the teal shading. The grey shading indicates the 1-σ confidence region of the sample of chunk measurements.

In the text

Fig. B.1. Summary of the method used to separate the two stellar populations along the RGB of M4. Left panel: The photometric data in the CU, B, I vs. V plane. Our targets are highlighted as orange points. The two dashed lines mark the location of the 10th and 90th percentiles of the RGB stars for the left and the right line, respectively. Upper right panel: The verticalised RGB (ΔCU, B, I vs. V) of M4. The points are divided between 1G (red) and 2G (blue) according to the GMM reported in the bottom right panel. The green points, seen far from the main RGB, represent the binary population. Bottom right panel: The GMM model used to identify the two stellar populations. The two most prominent Gaussian curves are the ones for the 1G and 2G stars.

In the text

	Fig. C.1. Correction values for our RGB targets (left), divided by populations, and our rHB stars (right). The solid lines are the corresponding models from which the corrections are drawn.
In the text

	Fig. E.1. Difference between the photometric (RCMD) and the asteroseismic (RS) radii for our RGB (top) and rHB (bottom) stars. The black dashed line in the upper panel is the average dispersion of the entire sample, while the coloured lines represent the average value of each population with their 1σ intervals as the shaded areas.
In the text

Fig. F.1. Systematic shift of the average mass and radii dispersion for different groups of stars. The first row reports the non-shifted values, which are also reported as transparent points in each panel. Each row past the first one explores the effects of a shift in a given parameter of its 1σ interval, as labelled. The last two rows report the non-corrected mass and dispersion average values, and what happens when the Temperatures coming from the B − V bands are used. The dashed lines represent the mass derived from the EB results in Kaluzny et al. (2013, 0.85M⊙, first column); the mass of the 1G, the 2G, and the HB from; Tailo et al. (2019, the values are 0.85M⊙, 0.833M⊙, and 0.64M⊙ for the second to fourth columns, respectively); and the shaded areas mark their 1σ intervals.

In the text

	Fig. F.2. Shift in temperatures for our RGB targets when the estimate is carried out with the B − V or the V − I colours. The dashed lines are the average values for the two populations with their respective 1σ interval.
In the text