Asteroseismic study of subgiants and giants of the open cluster M67 using Kepler/K2: expanded sample and precise masses

ABSTRACT

Sparked by the asteroseismic space revolution, ensemble studies have been used to produce empirical relations linking observed seismic properties and fundamental stellar properties. Cluster stars are particularly valuable because they have the same metallicity, distance, and age, thus reducing scatter to reveal smoother relations. We present the first study of a cluster that spans the full evolutionary sequence from subgiants to core helium burning red giants using asteroseismology to characterize the stars in M67, including a yellow straggler. We use Kepler/K2 data to measure seismic surface gravity, examine the potential influence of core magnetic fields, derive an empirical expression for the seismic surface term, and determine the phase term |$\epsilon$| of the asymptotic relation for acoustic modes, extending its analysis to evolutionary states previously unexplored in detail. Additionally, we calibrate seismic scaling relations for stellar mass and radius, and quantify their systematic errors if surface term corrections are not applied to state-of-the-art stellar models. Our masses show that the Reimers mass-loss parameter cannot be larger than |$\eta \sim 0.23$| at the |$2\sigma$| level. We use isochrone models designed for M67 and compare their predictions with individual mode frequencies. We find that the seismic masses for subgiants and red giant branch stars align with the isochrone-predicted masses as per their luminosity and colour. However, our results are inconsistent with the mass of one of the stellar components of an eclipsing binary system near the cluster turnoff. We use traditional seismic |$\chi ^2$| fits to estimate a seismic cluster age of |$3.95 \pm \, 0.35\, \mathrm{Gyr}$|⁠.

1 INTRODUCTION

Open star clusters are important laboratories for studying stellar astrophysics due to the shared properties among stellar members. The open cluster M67 is particularly important because of its solar-like age and composition, making it ideal for exploring the physics of stars that are only somewhat different to the Sun, which is the anchor point of most stellar models. This cluster is well-populated, providing a diverse and extensive sample of stars to study, which allows for thorough testing of stellar evolution theory (Gilliland et al. 1993).

Asteroseismology, the study of stellar oscillations, has firmly established itself as a powerful tool to study stars (Aerts 2021). The technique provides several advantages, such as the ability to quickly and accurately determine bulk properties, including surface gravity of large stellar samples across the Galaxy, which has provided a dramatic boost to the field of Galactic archaeology (Casagrande et al. 2016; Montalbán et al. 2021; Stokholm et al. 2023). Additionally, it has been suggested that asteroseismology can provide evidence of magnetic fields in the cores of red giants. Many observed red giants exhibit a low visibility of dipole oscillation modes compared to radial modes. This ‘mode suppression’ is theorized to result from the dissipation of dipole mode energy, caused by strong magnetic fields in the core. (Fuller et al. 2015; Stello et al. 2016b).

However, asteroseismology is often equated to just the application of seismic scaling relations (Brown et al. 1991; Kjeldsen & Bedding 1995). These relations, anchored to the Sun, rely on only two seismic observables (⁠|$\nu _{\mathrm{max}}$| and |$\Delta \nu$|⁠) to describe the global characteristics of the oscillations, providing estimates of stellar radius and mass, and, with the aid of models, stellar age. Despite their widespread use, the accuracy of these scaling relations, which are partly empirical in nature, remains not fully established (Huber et al. 2011; Miglio et al. 2012; Belkacem et al. 2013; Brogaard et al. 2018; Hekker 2020).

Beyond the scaling relations, the accuracy can be considerably improved, particularly of mass and age estimates, by measuring and performing detailed modelling of individual oscillation mode frequencies (Li et al. 2020; Miglio et al. 2021; Aguirre Børsen-Koch et al. 2022). Such ‘full’ asteroseismic analysis does, however, require corrections to the model frequencies because of improper prescriptions of the near-surface layers in current stellar models. As pointed out by Chaplin & Miglio (2013), these corrections should also be applied to the |$\Delta \nu$| scaling relation. Empirical methods exist to perform these corrections on a star-by-star basis (Ball & Gizon 2014), but unless very good seismic data are available, the corrections can introduce star-to-star biases. If these biases can be overcome, one interesting – and still unexplored – aspect of a full asteroseismic analysis is to investigate whether individual mode frequencies of both subgiants and red giants within a single stellar cluster can be consistently aligned with isochrone models that are matched to the classical observables.

M67 presents a unique opportunity to investigate the key aspects highlighted above. In particular, (1) To determine if seismic |$\log g$| aligns with state-of-the-art isochrones; (2) To assess to what degree mode suppression in this cluster follows the field population and whether clump stars are more or less susceptible to mode suppression than red giant branch stars; (3) To measure individual mode frequencies to obtain a more robust surface correction that evolves smoothly and map how the oscillation phase term |$\epsilon$| (of the asymptotic relation for acoustic modes) evolves during the poorly studied transition between subgiants and red giants; (4) To use individual frequencies to calibrate the scaling relations used to determine mass and radius; and (5) To employ individual frequencies for isochrone fitting.

In this study, we tackle all these aspects. We build upon the asteroseismology-based work on M67 presented by Stello et al. (2016c) and its follow-up study modelling three stars by Li et al. (2024). Both studies were based on single-campaign light curves from K2 (Howell et al. 2014) and focused solely on red giants. We expand on their research by including data from two additional K2 campaigns. The resulting improvement in signal-to-noise ratio led to the discovery of additional cluster members showing solar-like oscillations down to the subgiant branch and allowed us to identify individual frequencies for a large fraction of the seismic sample. We use the results of Reyes et al. (2024, R24), particularly their M67-specific isochrone models, distance modulus, and membership sample. This incorporates four additional oscillating red giants into the cluster membership over the membership sample used by Stello et al. (2016c) (The Geller, Latham & Mathieu 2015 sample). We focus on integrating seismic data from our extended sample with the isochrone models from R24 to derive highly precise stellar properties. In the following sections, we detail our methodology, present our findings, and discuss their implications.

2 SEISMIC ANALYSIS

2.1 Seismic data

We performed seismic analysis using time series photometry from K2 campaigns 5, 16, and 18, although not all stars have data in all three campaigns. Most stars were observed in long-cadence mode (⁠|${\sim} 30\, \mathrm{min}$|⁠) while the least evolved stars were observed in short cadence (⁠|${\sim} 1\, \mathrm{min}$|⁠). Long- and short-cadence light curves were created as described in Vanderburg & Johnson (2014), Stello et al. (2016c), and Vanderburg et al. (2016), except for four stars: EPIC IDs 211443624, 211414351, 211384259, and 211406144, which we downloaded from MAST (https://archive.stsci.edu/) as PDC-SAP flux corrected light curves.

The light curve preparation was performed as in Stello et al. (2016c) and included gap filling of up to three consecutive data points and high-pass filtering with a cut-off frequency of 3|$\,\mu$|Hz, except for stars with |$\nu _{\mathrm{max}}$||$< 8$| |$\mu$|Hz. We combined short cadence time series from different campaigns using statistical weights that take into account the varying data quality as a way to give lower weight to noisier parts of the time series. To calculate the weights, we used the rms point-to-point variability of the light curves in a 24-h rolling window. We then scaled the weights to be consistent with the noise measured in the amplitude spectrum, as described in Stello et al. (2006). Next, we calculated the weighted power spectrum as described by Frandsen et al. (1995) and Handberg (2013). For short-cadence data, the weighting scheme led to an average 12 per cent improvement in the signal-to-noise ratio in power. However, when we applied the same scheme to the long-cadence data, we did not see a significant improvement. We therefore decided to proceed without weighting the long-cadence light curves.

Upon analysing the oscillation spectra, we found that 51 stars showed solar-like oscillations, as indicated by a characteristic broad hump of excess power. As expected (Kjeldsen & Bedding 1995), the power excess is less pronounced in the highest |$\nu _{\mathrm{max}}$| cases, leading to a reduced signal-to-noise ratio. Despite this, we remain confident in our detections, as the power excess appears at frequencies consistent with predictions based on the stars’ positions in the colour–magnitude diagram.¹ Ten example stars, representative of the range of frequencies detected, are shown in Fig. 1. The power spectra of the seven M67 red clump stars are shown in Fig. F1. In Section 2.4, we outline our detection criteria for individual mode identification.

Among the 51 detections is EPIC 211 415 650 [S1072 from the Sanders (1977) catalogue, or WOCS 2008], a known evolved blue straggler, also known as a yellow giant, which is a member of a binary system (Mathieu & Latham 1986; Liu et al. 2008; Geller et al. 2015; Bertelli Motta et al. 2018). Evolved blue stragglers can be the product of stellar mergers or of mass-transfer episodes due to Roche lobe overflow of close binaries and do not follow the same standard evolutionary path as most stars in the cluster. EPIC 211415650, in particular, has been suggested to have formed through Kozai cycles (Perets & Fabrycky 2009) in a hierarchical triple system (Bertelli Motta et al. 2018). Therefore, we excluded EPIC 211 415 650 from the sample to be matched to our model isochrone, but we present its power spectrum and a complete seismic characterization in Appendix  A for future reference.

We used pysyd² (Chontos et al. 2021) to measure |$\nu _{\mathrm{max}}$| and |$\Delta \nu$| for our remaining sample of 50 M67 stars that show solar-like oscillations. In the following, we present results based on those |$\nu _{\mathrm{max}}$| and |$\Delta \nu$| measurements.

2.2 Temperatures and seismic log g

Fig. 2(a) shows Isochrone A from R24 and Gaia DR2 photometry (Gaia Collaboration 2018) of the stars in our seismic sample. The observed magnitudes have been converted to absolute magnitudes using a distance modulus of 9.614 calculated from Gaia DR3 (Gaia Collaboration 2023) parallaxes with zero-point corrections from Lindegren et al. (2021) and star-by-star corrections calculated from reddening dust maps (Planck Collaboration XIII 2016) (see R24).

To test whether the seismology agrees well with Isochrone A and therefore, with external observables of the cluster, we first calculated effective temperatures using the iterated InfraRed Flux Method (IRFM) and colour–|$T_{\mathrm{eff}}$| relations from Casagrande et al. (2021) using visible and near-infrared photometry from Gaia DR2 and 2MASS (Skrutskie et al. 2006). We adopted a metallicity of |$\mathrm{[Fe/H]}=0.05$| dex for all stars, to be consistent with Isochrone A. We then derived seismic |$\log g$| values using the expression |$g \simeq g_{\odot } (\nu _{\mathrm{max}}/\nu _{\mathrm{max},\odot }) \sqrt{(}T_{\mathrm{eff}}/T_{\mathrm{eff}, \odot })$| (Brown et al. 1991; Kjeldsen & Bedding 1995), with solar values of |$\log g_{\odot }$| = 4.4383 (Cox 2000) and |$\nu _{\mathrm{max}, \odot }$| = 3090 |$\mu$|Hz, and a first-iteration |$T_{\mathrm{eff}}$| from APOGEE DR17. In cases where the APOGEE temperature was not available, we used the GALAH DR3 (Buder et al. 2021) or Gaia DR2 temperatures in the first iteration. We obtained stable temperatures after three iterations and minimized their difference to Isochrone A by applying a |$-36.4 \mathrm{K}$| offset determined through least-squares optimization. It is not unusual to find systematic offsets of such magnitude when estimating effective temperatures. For example, Handberg et al. (2017) reported offsets of |$+34, +15, \mathrm{and} -45\, \mathrm{K}$| between their IRFM temperatures of the open cluster NGC 6819 and the averages of three other studies.

We compare our seismic |$\log g$| and IRFM |$T_{\mathrm{eff}}$| with Isochrone A in the Kiel diagram in Fig. 2(b), which shows a remarkable good agreement. We indicate non-single stars and stars with Gaia renormalized unit weight error (RUWE) exceeding 1.2 (Khan et al. 2019) with an open symbol. Non-single stars also include those identified by Geller et al. (2015) as binary, likely binary, stragglers, or possibly contaminated by a nearby source, and stars from the Gaia DR3 non-single star catalogue Gaia Collaboration (2023). These stars are not necessarily expected to closely follow the isochrone sequence, although most of them do. Additionally we show the primary component of the eclipsing binary WOCS1028 (Sandquist et al. 2021). The temperatures and |$\log g$|⁠, as well as |$\nu _{\mathrm{max}}$| and |$\Delta \nu$| for the full M67 sample, are presented in Table B1.

2.3 Mode suppression

Performing our initial seismic analysis, we found several stars whose |$\ell$| = 1 and/or |$\ell$| = 2 modes were undetectable or barely detectable. Therefore, we decided to measure their visibilities following Stello et al. (2016b), who found that the prevalence of mode suppression, thought to be due to interior magnetic fields (Fuller et al. 2015), is a strong function of stellar mass. Because our stars are in a cluster, we can unambiguously identify which of them belong to the helium-core burning clump; a task that is challenging for field stars with suppressed dipole modes. This difficulty led Stello et al. (2016b) to exclude red clump stars and focus solely on red giant branch stars observed by Kepler in long-cadence mode. In contrast, we measured visibility in all stars within our seismic sample with |$8.5<\nu _{\mathrm{max}}< 600\, \mu \mathrm{Hz}$| as this frequency interval offered both the resolution and signal-to-noise ratio necessary for a reliable measurement of the phase term |$\epsilon$| and therefore, a meaningful analysis of the visibilities. See Figs C1 and C2 for examples of |$\ell$| = 1 suppressed and non-suppressed red-giant-branch stars in M67.

Fig. 3(a) shows |$\ell$| = 1 visibilities from our M67 sample, where 11 of 29 red-giant-branch stars, and 1 of 7 red clump stars, are below the fiducial line from Stello et al. (2016b) that was found to separate the population of dipole suppressed red giants from the normal visibility stars. Fig. 3(b) shows the |$\ell$| = 2 visibilities, with an open symbol for stars with normal |$\ell$| = 1 visibilities and a filled symbol for |$\ell$| = 1 suppressed stars. We found that |$\ell$| = 1 suppressed red giant branch stars are more likely to also have lower |$\ell$| = 2 visibilities than normal red giant branch stars, in agreement with Stello et al. (2016a). We do not observe the same in red clump stars, but with only one dipole-suppressed star, the sample is too small to draw any conclusion in this regard. The low dipole suppression rate of |$1/7$| that we observe in red clump stars seems to corroborate the ‘most likely scenario’ proposed by Cantiello, Fuller & Bildsten (2016) for magnetic mode suppression in clump stars: that dipole suppression is rare in stars less massive than 2.1 M|$_\odot$| because strong core magnetic fields tend not to survive helium flashes.

To be able to compare the |$\ell$| = 1 suppression rate in M67 directly with that of the Kepler field stars measured by Stello et al. (2016b), we selected only red-giant branch stars with |$70<\nu _{\mathrm{max}}< 250 \,\mu \mathrm{Hz}$|⁠. From this selection, 6 of the 11 stars were |$\ell$| = 1 suppressed. This results in a fraction of |$0.55 \pm 0.15$|⁠, with uncertainty following binomial statistics: |$\sigma =\sqrt{p(1-p)/N}$|⁠, where p is the fraction of dipole-suppressed stars and N is the size of the sample (see Fig. 3c). Our results tend to indicate that the suppression rate in M67 exceeds that found in comparable field stars. This remains true if we remove the two suppressed stars closest to the fiducial line in Fig. 3(a), (⁠|$0.44 \pm 0.14$|⁠), or if we count them as non-suppressed stars (⁠|$0.36 \pm 0.15$|⁠). It would be interesting to study other Galactic clusters to find whether suppression rates are larger in cluster stars and what the implication might be for the magnetic activity of stellar cores in a cluster environment.

2.4 Individual radial mode frequencies

We expand the seismic analysis beyond just |$\nu _{\mathrm{max}}$| and |$\Delta \nu$| for both subgiant and red giant branch stars. We identified individual radial mode frequencies using the famed pipeline, which is based on Bayesian multimodal fitting and significance tests obtained from a model comparison process (Corsaro, McKeever & Kuszlewicz 2020). We note that famed is not optimized for our M67 K2 data, including aspects such as highly evolved stars, complex spectral windows, relatively short light curves, K2-specific noise characteristics and, for some stars, low signal-to-noise ratio. As a result, the Background tool facility diamonds + background (Corsaro & De Ridder 2014) meant to be used with famed³ often failed to converge, so we used the background models obtained by pysyd instead. The procedure used by pysyd to subtract the background of the power spectrum involves testing background models composed of one, two, or three Harvey components (Harvey 1985), with a white noise component either fixed or variable. It then computes the Bayesian Information Criterion (BIC; Schwarz 1978) for each model using the smoothed power spectrum and selects the best-fitting model based on the BIC value. On occasions famed mislabelled the degree of the modes or its process got aborted before all frequencies were scanned. In those cases, we corrected the degree or identified missing modes visually with the help of the echelle diagram. For the modes identified visually, we adopted the uncertainties given by famed to modes of similar signal strength from the same star or from one of the other stars in a similar phase of evolution. To ensure sample quality, we removed peaks with signal-to-noise ratio |$\mathrm{SNR}< 3.0$|⁠. However, we did not do this for stars with low overall signal-to-noise ratio; these are stars that we want to keep in the sample because their radial ridge was clear and their phase term |$\epsilon$| followed the cluster sequence, as we describe in Section 2.5. From our seismic sample of 50 stars, we excluded four highly evolved giants due to insufficient frequency resolution, one subgiant due to failure to identify oscillation modes confidently, and the seven helium-burning stars because modelling this evolutionary phase is beyond the scope of this work. In Table B2, we present the radial mode frequencies of the resulting sample of 38 stars. We used an asterisk (⁠|$^{*}$|⁠) to mark stars with a low overall SNR.

2.5 Surface correction

We followed the works by Kjeldsen, Bedding & Christensen-Dalsgaard (2008) and Ball & Gizon (2014), where the authors investigate the systematic difference, known as the surface term, between observed solar oscillation frequencies and those derived from a solar model. Instead of stellar models calibrated to the Sun, we use Isochrone A, which was calibrated to M67. And instead of a focus on how the surface term varies across radial mode orders of the Sun, we looked at the correction at |$\nu _{\mathrm{max}}$| and how that varies across many stars; those from Section 2.4.

The surface term is known to be caused by poorly modelled physical processes in stellar atmospheres (Christensen-Dalsgaard, Dappen & Lebreton 1988; Christensen-Dalsgaard et al. 1996; Houdek & Dupret 2015). Because real and model atmospheres change smoothly along the isochrone, we can expect a corresponding smooth change in the surface term (Li et al. 2023). This suggests that it should be possible to parametrize the M67 surface term, similar to the approach taken for the Sun by Kjeldsen et al. (2008) and Ball & Gizon (2014). To accomplish this, we obtained atmosphere-inclusive structure profiles from models along Isochrone A, with initial masses ranging between 1.30 and 1.38 M|$_\odot$|⁠, which span the masses from the main-sequence turn-off to the tip of the red giant branch.

We calculated adiabatic frequencies from the structure profiles using the oscillation code gyre (Townsend & Teitler 2013) version 6.0.1, scanning for radial modes up to the acoustic cut-off frequency |$\nu _{\mathrm{ac}} = \nu _{\mathrm{ac},\odot } \cdot (g/g_{\odot })/\sqrt{T_{\mathrm{eff}}/T_{\mathrm{eff},\odot }}$|⁠, taking |$\nu _{\mathrm{ac},\odot }= 5 \,\mathrm{mHz}$|⁠. For each star in Table B2, we selected the model that aligned best with the lowest order observed radial mode. Due to the low number of excited orders in some stars, we performed this selection manually noting that the lowest order mode is already likely affected by the surface effect in stars with few excited modes. We corrected the model frequencies for the surface effect by fitting the individual mode frequencies with the cubic expression described by equation (3) from Ball & Gizon (2014): |$\delta \nu _{\mathrm{corr}}(\nu) = a_{3}(\nu /\nu _{\mathrm{ac}})^{3}/I(\nu)$|⁠, where |$a_{3}$| is calculated as in their equations (8) and (9); and with the inverse-cubic expression |$\delta \nu _{\mathrm{corr}}(\nu) = [a_{-1}(\nu /\nu _{\mathrm{ac}})^{-1} + a_{3}(\nu /\nu _{\mathrm{ac}})^{3}]/I(\nu)$|⁠, where |$a_{-1}, a_{3}$| are calculated as in their equations (10) and (11). In both cases, |$\nu$| is the mode frequency, |$\nu _{\mathrm{ac}}$| is the model’s acoustic cut-off frequency, and |$I(\nu)$| is the mode inertia.

Fig. 4 shows in black symbols the resulting cubic surface corrections evaluated at |$\nu _{\mathrm{max}}$|⁠, as a function of |$\nu _{\mathrm{max}}$| across the M67 sample. A power-law fit to the cubic trend results in

as shown with a dotted line. The data point for the Sun was not included in the fit, but turns out to match our functional form when extrapolated. The resulting fit given by equation (1) was almost identical whether we used the cubic or the inverse-cubic expressions (grey-lime circles) from Ball & Gizon (2014), but the scatter around the fit was somewhat smaller in the cubic case.

Different model physics and in particular, different atmosphere choices can lead to drastically different values for the surface effect at |$\nu _{\mathrm{max}}$| (Christensen-Dalsgaard et al. 1996; Li et al. 2023). Therefore we note that our equation (1) is expected to be valid for models such as the ones described in detail in R24, which use the generalized Hopf-functions for atmosphere by Trampedach et al. (2014) as implemented by Ball (2021), with a modified radiative gradient that follows this T|$(\tau)$| relation.

In the phase-term diagrams of Fig. 5, the solid black curves show |$\Delta \nu$| and |$\epsilon$| of the models along Isochrone A, where we calculated (⁠|$\Delta \nu$|⁠, |$\epsilon$|⁠) as in White et al. (2011). We obtained the surface-corrected version of the isochrone models (dashed curves) by first deriving the surface term at |$\nu _{\mathrm{max}}$| from equation (1). We then substitute this value into the expression for |$\delta \nu _{\mathrm{corr}}(\nu)$| at |$\nu =\nu _{\mathrm{max}}$| to obtain the cubic coefficient |$a_{3}$|⁠. This coefficient allows us to calculate corrections for all radial modes in a given model. We repeated the procedure for all models along the isochrone to obtain the corrected (⁠|$\Delta \nu$|⁠, |$\epsilon$|⁠). The correction method had the effect of shifting |$\epsilon$| toward higher values relative to the uncorrected isochrone.

In Fig. 5(a) the black circles represent cubic-expression surface-corrected models found to be the best aligned to each of the stars, as described above. These individually corrected models show |$\epsilon$| values largely consistent with the surface-corrected isochrone and their scatter mirroring their scatter pattern seen around the fitted power law in Fig. 4. The dashed curve reflects the smoothing effect of the surface correction procedure implemented using the power-law fit. The extremely tight ‘evolution’ of |$\epsilon$| seen here has never been mapped out before due to the previous lack of very similar stars among field star samples and more generally the lack of stars with measured |$\epsilon$| in the transition between main sequence/early subgiant (Lund et al. 2017) and red giant branch stars (Corsaro et al. 2012), as nicely spanned by this M67 sample.

Fig. 5(b) follows the same concept, but with the grey-lime circles representing the inverse-cubic form of the surface correction to the same best-aligned models. Comparing Figs 5(a) and (b), it becomes clear that the cubic expression led to a smoother progression of the phase term |$\epsilon$| with |$\Delta \nu$| when compared to the inverse-cubic expression. The extra scatter observed with the inverse-cubic expression could be reflecting a natural scatter around the power-law fit if in fact the inverse-cubic corrections are more accurate. However, if the surface correction using the cubic expression is more accurate, this extra scatter could suggest that the inverse-cubic expression slightly overfits the data when, for example, one or more radial modes deviate from their expected frequency due to actual structural discontinuities or if the error in the peakbagged frequencies has been underestimated. In Fig. D1 we present the four stars whose surface-corrected |$\epsilon$| show the largest differences between using the cubic versus inverse-cubic expressions. Given that a smoother progression of the phase-term aligns better with the expectation from a uniform stellar sample such as ours, going forward we use the cubic fit.

Next, we calculated |$\Delta \nu$| and |$\epsilon$| using the observed radial mode frequencies (Table B2) and their uncertainties as weights, and we also used the observed radial modes along with pysyd-derived |$\Delta \nu$| to obtain pysyd-derived |$\epsilon$| values. These results are shown in Figs 5c and d, respectively. We observe a greater scatter in pysyd based |$\epsilon$| compared to the scatter in Figs 5a, b, and c, which we attribute to the effects of strong radial modes dominating over weaker ones, and to the effects of mixed modes, both of which influence |$\Delta \nu$| measurements made by autocorrelation-based pipelines like pysyd. These effects are amplified in |$\epsilon$| due to its high sensitivity to even small changes in |$\Delta \nu$|⁠. In Figs 5c and d, we observe a trend toward larger |$\epsilon$| compared to the model-derived |$\Delta \nu$| in Figs 5a and b, particularly for stars on the red giant branch. This discrepancy may result from the curvature of the radial ridge known to be more pronounced on red giant branch stars (see fig. 4 of Mosser et al. 2013). The curvature’s influence is more apparent in |$\Delta \nu$| from models, where there are many more orders available. This suggests that the weight function used to obtain |$\Delta \nu$| from models might not be narrow enough.

Fig. 6 illustrates the fractional differences between the observed |$\Delta \nu$| and |$\nu _{\mathrm{max}}$| values obtained from pysyd and those obtained through model matching. The figure suggests that observational features can influence |$\Delta \nu$| and/or |$\nu _{\mathrm{max}}$|⁠, but not necessarily both at the same time or to the same extent. If we assume that the model values are accurate, the figure shows that pysyd uncertainties are sometimes underestimated, not always reflecting the actual measurement errors. We also see that the fractional difference between the observed and model values tends to be higher for |$\nu _{\mathrm{max}}$| than for |$\Delta \nu$|⁠, suggesting that |$\nu _{\mathrm{max}}$| estimations based on the Gaussian fitting to the smoothed spectrum tend to be less robust than |$\Delta \nu$| estimation when the power distribution in the spectra is uneven, which can happen when the light curves are short.

The observed offset in |$\epsilon$| between the uncorrected isochrone and the surface-corrected isochrone shown in Fig. 5, and more explicitly in Fig. 7(a), confirms previous observations made by White et al. (2011) from stellar tracks and solar-metallicity field stars. Furthermore, our equation (1) provides a recipe for estimating such an offset in models of other stars, at least those with masses and metallicity comparable to M67, and provided that models similar to those of R24 are used. This was possible due to the size and ample |$\nu _{\mathrm{max}}$| range of the sample and the fact that the stars differ only (gradually) by small amounts in mass.

Perhaps more significant is the effect of surface corrections on |$\Delta \nu$|⁠. To quantify its impact on our isochrone models, we calculate the fractional difference between |$\Delta \nu$| from the uncorrected and surface-corrected models. This effect is particularly pronounced at low |$\nu _{\mathrm{max}}$|⁠. Given that in the widely used solar–scaling relation for mass (equation 5) |$\Delta \nu$| is raised to the power of −4, the fractional error shown in Fig. 7(b) becomes considerable, even for higher |$\nu _{\mathrm{max}}$|⁠. This underscores the importance of accurately accounting for surface corrections when interpreting modelled |$\Delta \nu$|⁠, particularly in the context of the mass scaling relations, which we examine in detail in the following section.

2.6 Calibrating mass and radius scaling relations

The asteroseismic scaling relations for |$\nu _{\mathrm{max}} \propto (M/R^{2})/\sqrt{T_{\mathrm{eff}}}$| and |$\Delta \nu \propto \sqrt{M/R^3}$|⁠) offer a simple way to determine stellar mass and radius (Stello et al. 2008; Kallinger et al. 2010). However, the relations are known to be only approximations (Stello et al. 2009; White et al. 2011; Miglio et al. 2013) and therefore corrections are needed to improve their accuracy (Sharma et al. 2016; Rodrigues et al. 2017; Kallinger et al. 2018; Hekker 2020). Our cluster stars provide a means to benchmark and potentially correct these relations, particularly those related to mass, which has been called into question in the past. For example, Rodrigues et al. (2017) suggested that we may have been overestimating scaling-based seismic masses, and Kallinger et al. (2018) said that scaling relations could not be applied directly to red clump stars, in agreement with Miglio et al. (2012) who had previously stated that due to the different temperature profiles between the red giant branch and the red clump stars, a relative correction has to be considered between both groups if we intend to use the |$\Delta \nu$| scaling to estimate the stellar mean density.

Following past cluster studies (Miglio et al. 2012, 2016; Handberg et al. 2017; Howell et al. 2022), we investigate scaling-based stellar masses using four combinations of seismic (⁠|$\Delta \nu$| and |$\nu _{\mathrm{max}}$|⁠) and classic (L and |$T_{\mathrm{eff}}$|⁠) observables, following the notation of Sharma et al. (2016):

where |$f_{\Delta \nu }$| and |$f_{\nu _{\mathrm{max}}}$| are correction factors. For instance, |$f_{\Delta \nu }$| has been calculated using model grids (Sharma et al. 2016; Rodrigues et al. 2017; Serenelli et al. 2017), and an expression for |$f_{\nu _{\mathrm{max}}}$| based on mean molecular weights has been suggested (Viani et al. 2017).

Similarly, radii will in the following be calculated with the seismic scaling relation:

except for the brightest star, EPIC 211376143, for which we adopt a radius of |$73.8\, \mathrm{R}_{\odot }$|⁠, as predicted by Isochrone A for the star’s |$\nu _{\mathrm{max}}$|⁠, because its |$\Delta \nu$| is unavailable due to insufficient frequency resolution of the data.

Because Isochrone A is fitted to Gaia DR2 magnitudes, we calculate the luminosities of our M67 stars using the same photometric system and the formula |$-2.5 \log L_{*}/L_{\odot } = G_{\mathrm{mag}} + BC_{\mathrm{ G}}-M_{\mathrm{bol},\odot }$|⁠, where |$G_{\mathrm{mag}}$| is the absolute G magnitude, and |$BC_{\mathrm{ G}}$| is the bolometric correction. Bolometric corrections depend on |$T_{\mathrm{eff}}$|⁠, |$\log g$|⁠, and [Fe/H], and we obtain them through interpolation to tables from the MIST project,⁴ which are based on the grid of stellar atmospheres and synthetic spectra described by Choi et al. (2016). |$M_{\mathrm{bol},\odot }$| is the absolute bolometric magnitude of the Sun that we take as 4.74 |$\pm$| 0.01 mag following Mamajek et al. (2015). To derive absolute G magnitudes, we use the distance modulus (⁠|$9.614 \pm 0.049$|⁠) and the differential extinction corrections derived in R24. To address the magnitude-dependent trend in DR2 G magnitudes seen by Evans et al. (2018) and Casagrande & VandenBerg (2018), we set the |$G_{\mathrm{mag}}$| uncertainties to 0.01 magnitudes, about 10 times the photometric precision estimated for DR2 sources at the brighter end (⁠|$G< 13$|⁠). Solar reference values are taken as |$\Delta \nu _{\odot }= 135.1 \,\mu \mathrm{Hz}$|⁠, |$\nu _{\mathrm{max}, \odot } = 3090\, \mu \mathrm{Hz}$|⁠, and |$T_{\mathrm{eff},\odot } = 5777 \,\mathrm{K}$| (Huber et al. 2011).

We calculated the masses and radii of the presumed single stars in our M67 sample, where the input |$\Delta \nu$| and |$\nu _{\mathrm{max}}$| was from our pysyd measurements, and |$T_{\mathrm{eff}}$| from Section 2.2 (Table B1), using equations (2) to 6. First, we assumed |$f_{\Delta \nu } = f_{\nu _{\mathrm{max}}}=1$|⁠. In Fig. 8, we compared these scaling masses and radii, represented by individual symbols, to the ‘true’ model masses and radii along Isochrone A (solid black curve). Here, each panel represents one of the four mass scaling relations. We found that the stars only align with the isochrone for the equation (2)-mass (the one without |$\Delta \nu$|⁠). In contrast, equations (3) to (5) show strong mass discrepancies that seem proportional to the power of |$\Delta \nu$| in each equation, suggesting a non-unity |$f_{\Delta \nu }$| factor is needed. To accurately address the mass discrepancies, which clearly vary in magnitude along the isochrone, it is essential to incorporate a |$f_{\Delta \nu }$| factor that depends on stellar parameters.

One common approach to determining |$f_{\Delta \nu }$| is by interpolating within pre-computed model grids, where |$f_{\Delta \nu }$| is defined as the ratio between |$\Delta \nu$| (measured from radial mode frequencies) and the square root of the mean stellar density. Examples include the grids by Rodrigues et al. (2017) and Stello & Sharma (2022). However, these grids are typically calibrated to the Sun and often do not account for surface effects. In Fig. 7(b), we previously demonstrated the impact of the surface term on |$\Delta \nu$|⁠. Fig. 9 demonstrates how ignoring the surface term will affect the mass estimates. The figure shows the scaling-based masses calculated using equations (2) to (5), similar to Fig. 8, but now incorporating |$f_{\Delta \nu }$|  = |$\Delta \nu _{\mathrm{obs}}/\sqrt{\rho _{\mathrm{model}}}$| (solar-scaled). This approach is similar to using the aforementioned grids, but with the added benefit that our isochrone is calibrated to M67. For each star, we determined the model density |$\rho _{\mathrm{model}}$| by interpolating the observed |$\Delta \nu$| (⁠|$\Delta \nu _{\mathrm{obs}}$|⁠, derived from pysyd and listed in Table B1) over the modelled |$\Delta \nu$| (⁠|$\Delta \nu _{\mathrm{model}}$|⁠), derived from radial frequencies along Isochrone A, prior to applying surface corrections. This approach represents an improvement over the case where |$f_{\Delta \nu }$| = 1, as the masses are now closer to agreeing amongst the four equations, and with the isochrone-predicted ‘truth’ values. However, significant discrepancies remain, in particular regarding the most commonly used equation (5). Fig. 9(d) shows mass discrepancies between |$\sim$|3 per cent and 20 per cent compared to the true isochrone masses. This agrees with the findings by Li et al. (2023). It highlights the need for a refined |$f_{\Delta \nu }$| approach, one that accounts for both the deviation of |$\Delta \nu$| from the square root of stellar density and the influence of the surface effect on |$\Delta \nu$|⁠.

We achieved this by applying |$f_{\Delta \nu }$| = |$\Delta \nu _{\mathrm{obs}}/\sqrt{\rho _{\mathrm{surf}}}$| (solar-scaled), where this time the density |$\rho _{\mathrm{surf}}$| is determined by interpolating the observed |$\Delta \nu$| over |$\Delta \nu _{\mathrm{surf}}$| from radial frequencies along Isochrone A after surface corrections. Fig. 10 shows that when using this more thorough correction, we find a very good alignment amongst the red-giant-branch scaling masses from the four different equations.⁵

To quantify and summarize the improvement of these corrections, we first show, in Fig. 11(a) (filled symbols), the average red giant branch mass across the black filled circles from Fig. 8 for each mass equation before applying any corrections. We then show the results after applying the square-root-density correction that neglects the surface term in Fig. 11(b). Comparing this with Fig. 11(c), which shows the averages after our last correction (square-root-density + surface term), clearly demonstrates that red giant branch scaling-based masses from each equation, on average, now agree with one another within uncertainties. Incidentally, they also align with Isochrone A’s mean red-giant-branch mass (solid black line), although that was not necessarily the intent of this exercise. However, even after this correction, in Fig. 11(c) the red clump masses (squares) still disagree amongst the different equations. Moreover, three of the averaged masses are so high that, if real, the stars would have evolved beyond the red clump phase by the age of this cluster. The exception this time is equation (3), the only one without |$\nu _{\mathrm{max}}$|⁠. We also see that the red clump mass disagreement seems to scale with the power of |$\nu _{\mathrm{max}}$|⁠. This suggests that this discrepancy can be removed using a |$\nu _{\mathrm{max}}$| scaling factor. As previously stated, the need for a red clump specific |$\nu _{\mathrm{max}}$| correction factor has been suggested before (Miglio et al. 2012; Viani et al. 2017; Kallinger et al. 2018), and in our data, a factor of |$f_{\nu _{\mathrm{max}}}=1.0477$| aligns all mean red clump masses to a consistent value of |$\overline{\mathrm{M}}_{\mathrm{RC}} = 1.387 \pm 0.044 \,\mathrm{M}_{\odot }$|⁠, as shown in Fig. 11(d). A summary of the scaling masses described here is presented in the table below Fig. 11.

The individual clump star masses that incorporate the |$f_{\nu _{\mathrm{max}}}$| correction are shown in Fig. 10 with grey squared symbols. Although the mass spread observed in red clump stars is still greater than that of the red giant branch sample in black circles, it is consistent with the mass spread observed in other lower SNR stars like the subgiants in our sample.

Fig. 11(d) shows that the averaged red-clump masses are compatible with the no-mass-loss scenario (solid grey lines in Fig. 11) in agreement with the no-mass-loss result – also from M67 asteroseismology – by Stello et al. (2016c). Specifically, we find the mass difference between red giant branch and red clump stars to be |$M_{\mathrm{RC}}-M_{\mathrm{RGB}} = 0.031 \pm 0.052\,$| M|$_\odot$|⁠, which should be compared to the no-mass-loss isochrone, which gives |$M_{\mathrm{RC,iso}}-M_{\mathrm{RGB,iso}} = 0.021\,$| M|$_\odot$|⁠. Hence, the maximum mass-loss consistent with our data within |$2\sigma$|⁠, |$(0.031-2\times 0.052)\,$| M|$_\odot$|⁠, would correspond to a Reimers |$\eta \sim 0.23$|⁠.

Interestingly, Pinsonneault et al. (2018) also found a negligible difference in mass between the RGB and RC in the metal-rich open cluster NGC 6791: |$0.02\pm 0.05\, \mathrm{M}_{\odot }$|⁠, using theoretically motivated corrections to the |$\Delta \nu$| scaling relations, in agreement with the results from Miglio et al. (2012), An et al. (2019), and Ash et al. (2025). These results suggest that mass-loss rates considered normal in the past (Reimers 1975; Carraro et al. 1996), may not be the norm in open clusters.

With the square-root density + surface term |$f_{\Delta {\nu }}$| corrections, and the red clump-specific |$f_{\nu _{\mathrm{max}}}$| correction factor that we found by unifying the mass agreement with all four mass scaling-based equations, we now want to confirm that those same corrections also improve the scaling-based radii (equation 6). Fig. 12(a) shows the fractional difference between the uncorrected scaling-based radii relative to the true isochrone radii, which for each star is defined as the radius of the model whose |$\Delta \nu _{\mathrm{surf}}$| matches the observed pysyd-|$\Delta \nu$|⁠. Fig. 12(b) shows that after applying our corrections, the agreement is indeed improved. Moreover, our best data, represented by the filled circles, show fractional differences remarkably close to zero. However, we see that some stars remain far from zero, despite having small uncertainties. This likely stems from underestimated |$\nu _{\mathrm{max}}$| and |$\Delta \nu$|  pysyd-uncertainties that occur due to the stochastic nature of the oscillations, as stated in Section 2.5 and also seen in Fig. 6.

To ensure that our mass results were not influenced by the models used to refine our calculations, we repeated our mass calculations for the best red giant branch data (filled black circles in Figs 8 to 10). We do so using our observed |$\nu _{\mathrm{max}}$|⁠, |$\Delta \nu$|⁠, |$T_{\mathrm{eff}}$|⁠, and L, but using the modifications to the mass equation (5) proposed by Li et al. (2022). This resulted in a mean red giant mass of |$\overline{\mathrm{M}}_{\mathrm{RGB}} = 1.347 \pm 0.072\, \mathrm{M}_{\odot }$| with their APOGEE-based fit, and |$\overline{\mathrm{M}}_{\mathrm{RGB}} = 1.383 \pm 0.076\, \mathrm{M}_{\odot }$| with their LAMOST fit. Our mean red giant mass estimate from the four equations |$\overline{\mathrm{M}}_{\mathrm{RGB}} = 1.356 \pm 0.027 \pm 0.008\, \mathrm{M}_{\odot }$| is consistent with both fits within the margin of error but is closer to their APOGEE fit, and has significantly lower uncertainties.

3 MODEL GRID FITTING

Having tested models from Isochrone A against M67 stars in the previous sections, we move on to perform model grid fitting using temperatures from Section 2.2, individual frequencies from Section 2.4, luminosities as described in Section 2.6, and models across an age range spanning 2 Gyr.

To construct the model grid, we evolved tracks with initial masses between 1.23 and 1.44 M|$_\odot$| with mass step |$10^{-4}$| M|$_\odot$| using the same model parameters as Isochrone A and saved atmosphere-inclusive structure profiles within a stellar age range of 3.0–5.0 Gyr. We calculated radial frequencies as in Section 2.5 for every model in the grid with |$\nu _{\mathrm{max}}$||$<$| 950 |$\mu$|Hz and an age gap between models of at most 2 million years, resulting in over two million individual models. As in Section 2.6, we used the power law derived in Section 2.5 to surface-correct the models’ radial modes using the cubic method, thereby obtaining a |$\Delta \nu _{\mathrm{surf}}$| for each model. For each presumed single star from Table B2, we created a subgrid with models whose |$\Delta \nu _{\mathrm{surf}}$| were within 1.5 times the average pysyd-|$\Delta \nu$| uncertainty of the observed |$\Delta \nu$| (⁠|$\sim$|1 per cent, Table B1).⁶ We then calculated a seismic |$\chi ^{2}$| from each star’s individual radial mode frequencies (Table B2), and models from its subgrid as follows:

where N is the number of observed radial modes for each star. We also derived a classical |$\chi ^{2}$| based on luminosity and effective temperature, assigning equal weights to both components:

Following Joyce & Chaboyer (2018) and to avoid introducing bias through arbitrary weighting choices, we combined seismic and classical ‘goodness of fit’ metrics into one |$\chi ^{2}$|⁠, as follows:

It is important to note that seismic and classical |$\chi ^2$| statistics are not independent, as both luminosities and temperatures are informed by the seismology. While we assign equal weights to these statistics, |$\chi ^{2}_{\mathrm{seismic}}$| imposes strong constraints, with only a narrow ‘slice’ of good-fit models along lines of equal density in the grid. Beyond this slice, the fit deteriorates rapidly as model frequencies deviate from observed values. In contrast, |$\chi ^{2}_{\mathrm{classical}}$| declines more gradually, introducing a gradient to the final |$\chi ^2$| statistic. This results in the best-fitting models being located within an elliptical region of the grid.

We selected a best-match result for each star based on the minimization of |$\chi ^2$|⁠. If the models in our grid accurately represent the seismic properties (stellar interiors) and classical observables of the cluster, the best-match models should align closely with the observations in the colour–magnitude diagram. This expectation is largely met in Fig. 13 (main panel), where Isochrone A is shown in black, and the full model grid is represented by dense grey dots, forming a broad grey band. While some individual matches are not perfect, importantly, there is no systematic offset (bias) where all the best-match models (red circles) are consistently brighter, fainter, redder, or bluer than the observed stars (white circles). This absence of directional bias suggests no overall density discrepancy between the matched models and the stars, a crucial observation for the discussion at the end of this section. This finding is further supported by the |$\chi ^2$| distribution along the grid for each star. The close-ups shown along the perimeter of Fig. 13 show each star’s subgrid, chosen based on the |$\Delta \nu _{\mathrm{surf}}$| proximity to the observed |$\Delta \nu$|⁠. The models from each subgrid are colour-coded according to their |$\chi ^2$| values and the diagonal stripes seen in most subgrids correspond to lines of equal density, with the darkest red stripe indicating models that are good matches to the observed mode frequencies and to |$T_{\mathrm{eff}}$| and L simultaneously.

The parameters of the best-fitting models are presented in Table B4. In Figs E1 to E4, we show the observed radial modes (small cyan circle) and best-fitting model (open orange circles) in echelle diagrams for each star. These figures include collapsed echelles to highlight the radial ridge positions.

According to |$\chi ^{2}$|⁠, the worse-fitted stars are those of EPIC IDs 211410817 (⁠|$\chi ^{2}$| = 1.1096), 211414300 (⁠|$\chi ^{2}$| = 1.1322), and 211413064 (⁠|$\chi ^{2}$| = 1.4938). Despite not being flagged as multiple star systems or having poor photometry, these stars deviate from the cluster sequence in the colour–magnitude diagram. This could indicate underestimated uncertainties, potential issues with their membership classification, unrecognized multiplicity, or other astrophysical phenomena affecting their photometry. Indeed, there is doubt about the membership classifications of 211 410 817 and 211413064. According to R24, these two stars were classified as cluster members by Geller et al. (2015), but could not be confirmed with Gaia DR3 astrometry because the radial velocities of 211 410 817 are not available, and the proper motions of 211 413 064 seem to contradict the membership status by Geller et al. (2015).

To get an overview of the fitting results for the sample as an ensemble, we show in Fig. 14 the ages and masses of the best-fitting model for each star, where the colour coding indicates the model’s |$\Delta \nu _{\mathrm{surf}}$|⁠. As expected, for a fixed age, the least evolved stars (dark blue) are of lower mass than the most evolved stars (dark red). The diagonal trend at fixed |$\Delta \nu _{\mathrm{surf}}$| (colour) reflects the intrinsic age–mass relation of the model grid, arising from age being the only parameter varied (because the other parameters were already fit to the cluster in R24). For clarity, we also show in this mass–age diagram all models in the grid, appearing as a continuous grey band. The histograms along both axis show the mass and age distributions of the |$\chi ^2$| fit.

Lastly, we estimated an average cluster age from the best-fitting models using the inverse of their combined |$\chi ^2$| as weights, represented by the circle size in Fig. 14. This resulted in a weighted average cluster age of |$3.95 \pm 0.35\, \mathrm{Gyr}$|⁠,⁷ where the uncertainty is derived from the standard deviation of the age values in the best-fitting models and is indicated by horizontal dashed lines around the average. This age aligns with that obtained via classical isochrone fitting to the cluster’s main sequence and red-giant branch stars by R24, though we obtain a larger uncertainty here. The fact that we obtain an age similar to R24 is likely because all models along the darkest red stripes in the subgrids of Fig. 13, which indicate low combined |$\chi ^2$|⁠, are also good seismic matches. There, younger models would appear on the left side of the stripe and older ones on the right). Thus, age is predominantly determined by |$\chi ^{2}_{\mathrm{classical}}$| along the stripe, closely approximating classical isochrone fitting. Our larger age uncertainty here may be attributed to our smaller sample size of 28 stars compared to the nearly 200 stars in the R24 sample, which also included main sequence and turn-off stars. Because we find the same age as R24, like them, we conclude our age estimate is in agreement with the estimates from VandenBerg & Stetson (2004) (⁠|$4.0 \pm 0.4\, \mathrm{Gyr}$| from |$BV$| photometry), Sarajedini, Dotter & Kirkpatrick (2009) (⁠|$3.5{\!-\!}4.0\, \mathrm{Gyr}$| from 2MASS photometry), Bellini et al. (2010) (⁠|$3.9 \pm 0.1\, \mathrm{Gyr}$| from the magnitude of the white-dwarf cooling sequence), Barnes et al. (2016) (⁠|$4.2 \pm 0.2\, \mathrm{Gyr}$| from gyrochronology), and Nguyen et al. (2022) (⁠|$3.9\, \mathrm{Gyr}$| from Gaia DR2 photometry using PARSEC isochrones v2.0).

The fact that our averaged age agrees with the age of Isochorne A means that Isochrone A agrees not only with the classical observables of all the cluster stars (as per R24) but also with the radial mode frequencies of the subgiants and red giants (this work). This has implications for the unresolved issue of the eclipsing binary WOCS 11028. Notably, Sandquist et al. (2021) reported precise masses for the eclipsing binary WOCS 11 028 near the cluster’s TAMS, but the primary component, WOCS 11028a, does not align with any existing isochrone for the cluster (Sandquist et al. 2021; Nguyen et al. 2022). There are discrepancies of at least 0.17 mag between the deconvoluted primary photometry and the predicted luminosity for the primary’s mass (⁠|$1.222 \pm 0.006$| M|$_\odot$|⁠), with the discrepancy with Isochrone A being 0.2 magnitudes. To further investigate this discrepancy, we applied our |$\chi ^2$| fitting approach as described above, but instead, we used models from an alternative isochrone made to fit the primary component and the rest of the cluster’s morphology in the colour–magnitude diagram (Isochrone EB from R24). In that case, the fit revealed a one-directional offset for all stars in our sample, where the best seismic fits on the red giant branch were at least 0.1 mag fainter than the observed values. This further confirms that the mass of WOCS 11028a cannot be reconciled with the stars in our seismic sample in a scenario where they all have evolved independently and were born from the same primordial molecular cloud. This therefore leaves us with two possibilities. Either the inconsistency stems from WOCS 11028a, or more concerning, it could suggest our current models may be flawed for masses (or evolutionary phases) close to the ones of WOCS 11028a. Given that (1) the cluster membership of WOCS 11 028 is undisputed, (2) there are indications the binary components have not interacted in the past (Sandquist et al. 2021), and (3) their fundamental parameters are fairly well established, our inability to align WOCS 11028a with the rest of the stars supports the need for further studies into the issue.

4 SUMMARY, CONCLUSIONS, AND FUTURE WORK

In this work, we used Kepler/K2 data to obtain seismic parameters of subgiants and red giants of the Open Cluster M67. Specifically:

• We obtained and evaluated |$\log g$| against the modelled isochrone from R24 in the Kiel diagram and found remarkably good agreement.

• We measured the ratio of cluster stars with possible internal magnetic fields in their cores using the mode suppression method and found that the suppression rate in M67 exceeds that found in comparable field red giant branch stars. By looking at the suppression rate of red clump stars, we provide evidence that magnetic cores do not normally survive the helium flash.

• We measured individual mode frequencies, which allowed us to develop an expression for correcting the surface term that changes smoothly with stellar evolution.

• For the first time, we mapped the evolution of the phase term |$\epsilon$| in the transition between main sequence turn-off and red giant branch stars.

• We successfully calibrated the masses and radii from scaling relations against our models and showed how reliable the masses derived from scaling relations truly are.

• We performed |$\chi ^2$| fitting over a grid of models using the radial modes of 28 stars, including giants and subgiants. This approach yielded the most precise subgiant and red giant branch star masses reported to date for the cluster. Importantly, our seismic fitting results are consistent with classical isochrone fitting methods. Still, our age estimates remain subject to potential age systematics stemming from our choices of input physics (see figs 42 and 43 in Lebreton, Goupil & Montalbán 2014). However, we note that the corrections applied in this work, whether for |$T_{\mathrm{eff}}$| offsets, the surface term, or |$f_{\nu _{\mathrm{max}}}$| and |$f_{\Delta \nu }$| correction factors agree with literature values. In this context it is also worth noting that our models agree with so-called small frequency separations, which probe the deep interiors of stars (Reyes 2025, in press).

The main avenues for future work include a more detailed investigation of the red clump stars, such as identifying their individual frequencies, conducting in-depth comparisons with models, and assessing whether the surface correction differs from that of the red giant branch stars. Furthermore, one could incorporate non-radial mixed modes in the analysis for all the stars. One could potentially obtain improvements by using SED fitting to derive effective temperatures and bolometric fluxes. Finally, the characterization of the yellow straggler (Appendix  A) would benefit from a comprehensive modelling approach that includes uncertainty estimations of the fundamental stellar properties.

SOFTWARE

This work uses the following software packages:

• famed version 0.0.1 Corsaro et al. (2020).

• mesa version 23.05.1 Paxton et al. (2011, 2013, 2015, 2018, 2019); Jermyn et al. (2023).

• gyre version 6.0.1 Townsend & Teitler (2013).

ACKNOWLEDGEMENTS

DS is supported by the Australian Research Council (DP190100666).

RDM acknowledges the support of NASA grants NNX15AW69G, NSSC18K0449, and NSSC19K0105.

The authors thank Joel Ong, Guy Davies, and Jim Fuller for useful comments on the manuscript.

This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/Gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/Gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This research includes computations using the computational cluster Katana supported by Research Technology Services at UNSW Sydney, https://doi.org/10.26190/669x-a286.

DATA AVAILABILITY

Isochrone A is available for download at DOI: 10.5281/zenodo.12616440.

The grid of models presented in this work can be found at DOI: 10.5281/zenodo.14927665.

Footnotes

REFERENCES

APPENDIX A: YELLOW STRAGGLER EPIC 211415650

We report the first seismic analysis of the yellow straggler star EPIC 211415650. This star had short cadence data from campaigns 16 and 18 only, the last data received from K2. We treated the data as described in the main text. We calculated all seismically derived parameters in the same way as for the main sample of this work, applying the described corrections to the scaling relations. A 1.84 M|$_\odot$|⁠, 1.48 gigayears old model calculated similarly to the models from Isochrone A is a good fit to the observed colour and magnitude of the yellow straggler and simultaneously to its observed |$\nu _{\mathrm{max}}$| and |$\Delta \nu$|⁠. See Fig. A1.

APPENDIX B: TABLES

APPENDIX C: VISIBILITY MEASUREMENTS

We present examples of visibility measurements presented in Fig. 3. Fig. C1 shows red-giant-branch stars of normal |$\ell$| = 1 visibility, and Fig. C2 shows examples of |$\ell$| = 1-suppressed red-giant-branch stars.

APPENDIX D: STARS WITH LARGE PHASE-TERM DIFERENCES: CUBIC VERSUS INVERSE-CUBIC

In Fig. D1, we present echelle diagrams for four stars, selected because they exhibit the largest differences in |$\epsilon$| between their cubic and inverse-cubic surface-corrected models within the sample. Each row corresponds to the same star, with the echelle diagrams constructed using either the cubic or inverse-cubic surface-corrected |$\Delta \nu$| and displaying the corresponding surface-corrected modelled frequencies, as indicated.

For star 211406541, neither the cubic nor the inverse-cubic correction appears optimal. This mainly arises because the number of observed modes is too low to describe the strong curvature in the radial ridge of the model with high confidence. As a result, the cubic fit (left) does not align well with the observed modes (in yellow), particularly around 30–40 |$\mu$|Hz. Although the inverse-cubic fit (right) does align with the observed modes, the reliability of the fit is questionable because the resulting |$\epsilon$| value falls far from the smooth cluster sequence in Fig. 5.

For star 211411922, it is also challenging to determine which is the better correction, this time because both models fit the observed frequencies reasonably well. Here, only three radial modes were identified and used for the fits, resulting in larger uncertainties and hence different |$\epsilon$| values.

For stars 211417812 and 211421683, the higher order radial modes in the inverse-cubic corrected models tend to fall in regions with no power beyond the fitted modes. For the cubic-corrected models, those modes coincide with potential low-amplitude modes (not highlighted in yellow). These interpretations are speculative, because the signals under consideration are statistically insignificant.

APPENDIX E: ECHELLE DIAGRAMS: STARS AND THEIR BEST-FITTING GRID MODELS

APPENDIX F: RED CLUMP STARS

Author notes