Stellar Evolution in Real Time II: R Hydrae and an Open-Source Grid of
$>3000$ Seismic TP-AGB Models Computed with MESA

The presence of TiO and VO molecular bands in the spectra of R Hya indicate that it is an O-rich Mira (Lockwood, 1969). Many other oxygen-bearing molecules have been observed in the circumstellar envelope of the star as well, as described more recently in Wallström et al. (2024).

The distance to R Hya was determined to be $125\pm 10$ pc by Haniff et al. (1995) based on period–luminosity (PL) relations. This agrees well with the result of Andriantsaralaza et al. (2022), who found $d=126_{-2}^{+3}$ pc based on a reanalysis of Gaia DR3 data. Andriantsaralaza et al. (2022) also derived a luminosity of $L=10300\pm 300\,L_{\odot}$ and effective temperature of $T_{\rm eff}=3100$ K using DUSTY models fitted to the spectral energy distribution (SED) of the star. Other authors found lower luminosity and temperature values based largely on broad-band near-infrared color data: De Beck et al. (2010) derived $L=7375\,L_{\odot}$ and $T_{\rm eff}=2128\pm 140$ K; Haniff et al. (1995) found $T_{\rm eff}=2680\pm 70$ K; Feast (1996) found $T_{\rm eff}=2830\pm 170$ K.

Given the proximity and large physical size of the star, its radius can be derived via interferometry. Feast (1996) published $R=401\pm 46\,R_{\odot}$, whereas Haniff et al. (1995) published two different model solutions corresponding to $R=442\pm 65\,R_{\odot}$ or $R=384\pm 54\,R_{\odot}$. We note that these measurements refer to the contemporary radius value and not to the maximum extent of the star when it was pulsating with longer periods.

The steady decrease of the pulsation period has been known for over a century, with regular observations going back to the 1840s. Zijlstra et al. (2002) successfully recovered the pulsation period for the second half of the 17th century and late 18th century from historical maps, and concluded that R Hya appeared to have a constant period of 495 d before the start of the period decline. Based on their TP models, Wood & Zarro (1981) proposed that R Hya is evolving from the maximum expansion and luminosity phase of the TP, and the change in $\dot{P}$ matched the change of slope in their model.

After the extensive decline, the period seemed to be stabilizing at 385 d in the 1990s. This raised the possibility of a mode switch due to envelope relaxation (Zijlstra et al., 2002), or an end to the radial decline (Fadeyev, 2023).

Since then, however, more than 20 years of additional photometry has become available. Figure 1 shows the recorded light variations of R Hya over 380 years, extending into 2023. In this plot, we combine the visual brightness estimates collected by the American Association of Variable Star Observers (AAVSO) with the measurements of the Digital Access to a Sky Century at Harvard project (DASCH¹¹1https://dasch.cfa.harvard.edu/) based on archival photographic plates (Grindlay et al., 2012). We scaled the average brightness and average pulsation amplitude of the DASCH data to match that of the visual estimates.

The AAVSO light curve immediately suggests that the pulsation amplitude has continued to decrease over the last 15–20 years. Our analysis of the new data set indicates that the period change of R Hya did not stop, but rather began to decline steeply again since 2000. Figure 2 shows the full span of observations, including the historic detections and period estimates compiled by Zijlstra et al. (2002).

We estimated the more recent periods in two ways for Figure 2. First, we cut the AAVSO and DASCH data into 2000 d long segments and calculated the pulsation frequencies of each segment with Period04 (Lenz & Breger, 2005). We fitted the pulsation frequencies along with two harmonics. For the AAVSO data, we used a sliding window with a step size of 1000 days. We also used the pulsation cycle lengths calculated from maximum and minimum timings published by AAVSO (2021) and by Karlsson (2013) in his extensive Mira O–C database²²2http://var.saaf.se/mirainfooc.php. The latter helpfully includes the 19th century observations collected in the seminal paper of pioneering female astronomer Annie Jump Cannon (Cannon & Pickering, 1909)³³3We encourage the reader to check this paper for bright long-period variables manually, because it is not indexed by SIMBAD. in digitized form. Unfortunately, the database comes without error values. The authors t hemselves note that in many cases they had no information on the uncertainties in the maxima timings, which could be as high as multiple days, depending on observers and observing conditions.

Based on the extensive period data, we find that, instead of stagnating in its period decline, the star is experiencing meandering, a well-known behavior of Mira stars in which the pulsation period changes quasi-periodically on a decadal to centennial scale. The combined effect of the meandering and a continuing evolutionary decrease led to an apparent but temporary stabilization of the pulsation period. The new observations extend the length of the period change to nearly 250 years and can be fitted with a linear rate of $\dot{P}=-0.493\pm 0.011$ d/yr.

The physical cause of meandering has not yet been established. One hypothesis suggests that meandering is caused by thermal relaxation oscillations, potentially connected to the TP (Templeton et al., 2005). However, current TP-AGB evolutionary models are unable to reproduce meandering; they are designed to calculate smooth changes in the stellar structure during oscillations caused by He-shell flashes that take place over few-hundred year timescales. Since the physics of meandering is not included in MESA, we can only fit the global period decline trend underlying the meandering with our models. This means that our theoretical curves may appear to be poor fits to some individual period measurements, and the goodness-of-fit calculations will be affected by meandering-based offsets that are not reflective of the true evolutionary trend.

Because of this, we tested whether separating the meandering from the TP-induced period changes resulted in better model fits. We did this by fitting a quadratic polynomial to the post-1845 period values. We fitted both the AAVSO light curve segment periods and the Karlsson (2013) data set: the differences were only noticeable for the very sparse late-19th century segment of the data. We elected to work with the AAVSO segments since we could determine period uncertainties for them, whereas the maximum timings did not include any uncertainty data. The best-fitting curve is the following: $0.00086819\,t^{2}-0.5839\,t+441.279$, where $t$ is defined in years relative to the epoch of JD2400000. With the presumed TP-induced period change removed, we can compare the meandering signal to that of other Miras. A frequency analysis reveals a dominant periodicity of 36.8 yr.

However, applying this fit involves making an assumption about the shape of the period decline. Instead, we adjust our fitting method to handle the elevated fit uncertainties coming from the presence of extra variations caused by physics that is not included in the models.

The pulsation amplitude of the star appears to be shrinking in parallel with its period. Over the last 130 years, enough photometric data have been collected to analyze the temporal evolution of the pulsation amplitude and average brightness of the star. We use the Fourier coefficients calculated in Section 2.2 and find that the Fourier amplitude of the pulsation frequency and the frequency itself are strongly correlated. As the frequency increases (period shrinks), the amplitude decreases as well.

We also examined the average brightness of the star. If R Hya is in a TP, its luminosity should be decreasing along with its radius. At first glance, the average brightness appears to be increasing in Figure 2, but that is only an artefact of using magnitude units. The visual brightness of the star changes by multiple orders of magnitude, and the bright phases get compressed by the logarithmic scaling. If we convert it into flux units, the downward change in pulsation-averaged brightness becomes apparent. We note that we did not apply bolometric corrections to the visual flux. The variation of Mira stars is much larger in visual than in bolometric units due to flux conversion into infrared near minimum light, therefore the average flux relative to solar we plot here is also lower than the luminosity of the star. Nevertheless, we observe a strong correlation with the average flux and the pulsation period.

Mira stars have cool atmospheres with a rich array of absorption features, including molecular lines. However, this makes determining their [Fe/H] exceedingly difficult. Here we can turn to clusters that contain Mira variables, and use the cluster metallicities as proxies. The period–[Fe/H] relation proposed by Feast & Whitelock (2000) suggests an [Fe/H] range of 0.1 to 0.65 dex, indicating super-solar metallicity. However, this is based on extrapolation, since the clusters used by Feast & Whitelock (2000) all have sub-solar metallicities.

The location and kinematic information of R Hya can further constrain its physical parameters. We calculated the galactocentric coordinates and velocities of R Hya using the Gaia DR3 data (Gaia Collaboration et al., 2023) and the python notebooks for the Gala python package (Price-Whelan, 2017; Price-Whelan et al., 2020). Since there is no radial velocity data for R Hya in DR3, we substitute the value of –15.9 km/s measured by Uttenthaler et al. (2011). The Galactic velocity components of the star, relative to the Local Standard of Rest (Schönrich et al., 2010), are $U_{\mathrm{LSR}}=-41.1\pm 1.0$ km/s, $V_{\mathrm{LSR}}=-9.24\pm 0.04$ km/s, and $W_{\mathrm{LSR}}=2.05\pm 0.04$ km/s, respectively. These low velocities place the star in the thin disk population on a Toomre diagram (e.g., Feltzing et al., 2003; Chen et al., 2021).

We find R Hya’s galactocentric distance and the vertical distance from the galactic plane to be $R=8041\pm 159$ pc and $z=114\pm 2$ pc, respectively. If we compare this to the observed metallicity distribution functions of the APOGEE survey, we find that the star falls into the group centered on solar metallicity, with the distribution extending from [Fe/H] = - 0.6 all the way up to +0.5 (lower panel of Fig. 5 in Hayden et al. 2015, $|z|<0.5$ kpc, $7<R<9$ kpc slice). This is consistent with the range we obtained from the period–[Fe/H] relation, and suggests that R Hya may indeed be super-solar in terms of metallicity.

R Hya is more accurately the R Hya system. It has one resolved companion with an estimated mass of 0.8 $M_{\odot}$ at a separation of 3500 au, which is too far away and too small to alter R Hya’s evolution (Mason et al., 2001; Anders et al., 2019). However, there are some indications of a much closer companion as well. High-resolution observations revealed bi-conical outflows around the star that is potentially shaped by a close-by companion (Decin et al., 2020). Homan et al. (2021) observed a disk around the star with an inner edge at a distance of 6 AU (or $\sim 3$ stellar radii). Based on the geometry and velocities of the disk, they estimate a current central mass of $\sim 2.5\,M_{\odot}$, which includes a proposed inner companion with a tentative lower mass limit 0.65 $M_{\odot}$.

The presence of the distant companion offers the possibility of measuring the metallicity of the system via that companion instead of R Hya itself. This star, Gaia DR3 6195030801634430336, is a K-type dwarf (Smak, 1964), and metallicities of (late) K- and M-type dwarfs could be estimated either from J–K photometry or from moderate-resolution spectroscopy (Johnson et al., 2012; Mann et al., 2013). Unfortunately, we did not find the appropriate observations. Nevertheless, Gaia DR3 6195030801634430336 offers an opportunity to verify the metallicities predicted by the AGB models for R Hya. We therefore urge observers to target this dwarf companion.

While the possibility put forth in Homan et al. (2021) is worthy of future consideration, theirs is currently the only claim of a close companion. We proceed treating R Hya using a single-star evolutionary assumption, which is also appropriate for well-separated, non-interacting binaries (e.g., Joyce & Chaboyer 2018). When we refer to R Hya henceforth, we are referring to the primary component of the wide binary.

The physical assumptions adopted in the evolutionary calculations include the Grevesse & Sauval (1998) opacities and solar mixture, an Eddington grey $T$–$\tau$ relation for the atmospheric boundary conditions (Eddington, 1930), the ‘Henyey’ mixing length scheme (Henyey et al., 1964) with a fixed value of $\alpha_{\text{MLT}}=1.931$ per the calibration to AGB stars performed in Cinquegrana & Joyce (2022a), the ‘pp_extras.net’, ‘co_burn.net’ and ‘approx21.net’ nuclear reaction network options, the Schwarzschild criterion for convective stability, and the ‘Reimers’ and ‘Blocker’ schemes for cool RGB winds and cool AGB winds, respectively (Reimers, 1975; Bloecker, 1995). Because we use the Schwarzchild criterion, we do not take thermohaline mixing or semi-convection into account. We do not consider rotation. We do not study the evolution of binaries or multiples but reiterate that single-star evolutionary calculations are acceptable for well-separated, non-interacting binaries like R Hya and its distant companion, Gaia DR3 6195030801634430336.

Though it is well known that the mass loss treatment in AGB models can significantly impact the number of pulses a star undergoes before exhaustion (e.g., Ventura & D’Antona 2005; Stancliffe & Jeffery 2007; Marigo et al. 2020), we studied only one value for the mass loss coefficient in the Blöcker wind scheme, $\eta_{\text{Bl\"{o}cker}}=0.01$. This value was chosen in part for consistency with the mass loss scheme and efficiency assumed in the Monash stellar evolution code (see Section 11) and to follow the “default” conventions assumed in, e.g., the MIST isochrones and stellar tracks (Choi et al., 2016). It is expected that lowering the mass loss efficiency would increase the total number of thermal pulses, and raising it would do the opposite. The anticipated effects of changing the mass loss prescription are less obvious. To the best of our knowledge, no systematic study of the impact of varying mass loss prescriptions in MESA has been conducted, though the literature would clearly benefit from such a study.

Likewise, while it was shown in Cinquegrana et al. (2022) that the AESOPUS opacities (Marigo & Aringer, 2009) are most appropriate for AGB stars, convergence difficulties and a desire for grid uniformity necessitated a reversion to the (better tested) Ferguson et al. (2005) low-temperature opacities, which were adopted throughout. Likewise to mitigate convergence difficulties, the models do not include convective overshoot. While this omission will produce slightly incorrect core masses, it does not affect our determinations of the best fit to R Hya, which are based on surface and near-surface parameters. However, we direct the reader to Cinquegrana & Joyce (2022a), Cinquegrana et al. (2022) and Cinquegrana et al. (2023) for best practices regarding the treatment of opacities, convective boundary mixing, and other physics for AGB models computed with MESA.

To keep computing time reasonable, the structural and time resolution coefficients, mesh_delta_coeff and time_delta_coeff, respectively, are not changed from their defaults.⁵⁵5 Except in the case of a handful of non-adiabatic calculations with GYRE, for which mesh_delta_coeff was halved, corresponding to a doubling of MESA’s structural resolution. The time resolution is not a concern, as the models are mathematically “well behaved” during the relevant regions and cubic splines are used for interpolation (see Section 6). However, improving the default structural resolution by a factor of 10, i.e., setting mesh_delta_coeff = 0.1, has been shown to produce changes in predicted $p$-modes of up to $0.4\,\mu{\rm Hz}$ on the main sequence (M. Joyce & Y. Li 2024, in prep). While the extent of this effect in supergiants is still being investigated, we introduce a 60 d buffer in period when we set the requirements for agreement between the observations and models (Table 2) to take this convergence-related modeling uncertainty into account.

The assumptions adopted for the adiabatic asteroseismic calculations did not deviate significantly from GYRE defaults. One notable difference is use of the ‘MAGNUS_GL6‘ shooting scheme, or sixth-order Gauss–Legendre Magnus method, in the numerical options. This is more reliable for finding modes in evolved stars as compared to the default diff_scheme option, which is second-order. The complete lists of parameter values for all MESA and GYRE controls are available in those files (i.e., inlists) on Github⁶⁶6 https://github.com/mjoyceGR/AGB_grid_visualizer and Zenodo⁷⁷7Zenodo link available upon publication.

An example of a seismic pulse spectrum produced with this configuration is shown in Figure 4. The periods, in units of days, are shown as a function of time during the TP-AGB phase. Each curve represents one of the 10 lowest-frequency (longest period) radial $p$-modes, $n_{p}=\{0,..,9\}$; $\ell,m=0$. An equivalent data set is available for every model in the grid.

Because the pulsation period measurements of R Hya constitute our primary observable, it is important to have an exact solution for the star’s radial oscillations at every time step in the evolutionary calculations. However, this is not an out-of-the-box capability of MESA itself; rather, MESA provides structural profiles that can be used as initial conditions for seismic calculations. The GYRE program generates a frequency spectrum characterizing a structural model’s response to perturbation, which can then be compared directly to asteroseismic observations.

It is computationally expensive to run an evolutionary model that associates a theoretical frequency spectrum to every time step for the entire duration of the TP-AGB, especially with the temporal resolution required to study thermal pulses in detail (5–10 years in the case of T UMi; Molnár et al. 2019). This issue was dealt with in Molnár et al. (2019), hereafter Paper I, by first manually isolating the regions of interest on evolutionary tracks, then re-running seismic calculations for the selected regions only. Among the many shortcomings of this technique are the lack of automation—making the computation of large grids intractable—and the generation intermediate structural data (profile) files that served as initial conditions for separate, corresponding GYRE runs.

We have improved upon the techniques of Paper I in three ways:

1. 1.

   We ran GYRE as part of the evolutionary calculations themselves rather than as a separate, post-processing step, thus decreasing the number of data output and writing operations (time⁸⁸8Due its greater parallelization, GYRE benefits more from multiple threads than MESA does. While low-order $p$-mode calculations are not particularly expensive, it is useful to allocate more threads to simulations calling GYRE on-the-fly. A model running on two threads that calls GYRE at every time step will run more slowly than the same model, running on the same number of threads, that does not. and data volume reduction);

2. 2.

   We relied on more sophisticated fitting and statistical techniques in lieu of greater time resolution (significant run time reduction); and

3. 3.

   We developed several interactive data visualization tools that led to improvement in the accuracy and scope of the results.⁹⁹9in Python, publicly available at the first author’s Github: https://github.com/mjoyceGR/AGB_grid_visualizer

Whereas Paper I treated MESA and GYRE calculations as independent, we used a technique called GYRE on-the-fly for this analysis. This method was adapted from a 2022 MESA Summer School group laboratory exercise authored by Earl Bellinger, available freely on Zenodo (earl_patrick_bellinger_2022_7118662)¹⁰¹⁰10Exercise instructions available here: https://earlbellinger.com/mesa-summer-school-2022/ Some of this functionality is also present in select MESA test_suite cases, available with the code, though documentation is more limited in those cases.

Calculations are broken into three evolutionary stages: ZAMS (zero-age main sequence) to TACHeB (terminal-age core helium burning), or phase 1, TACHeB to the onset of the TP-AGB phase (phase 2), and TP-AGB through termination (phase 3). Termination is determined by one of the following:

• -

  reaching the end of core helium exhaustion and ascent onto the WDCS (White Dwarf Cooling Sequence), required by the termination string ‘phase_WDCS’;

• -

  the condition min_timestep_limit, which is an indication of convergence failure;

• -

  the condition max_number_retries $=500$, which is an indication of convergence difficulty; or

• -

  run time exceeding 72 hours on 16 cores ($Y_{i}$-fixed grid) or 120 hours on 16 cores ($Y_{i}$-varied grid).

The latter two conditions are practical (not physical) choices that balance the goal of fitting R Hya specifically—for which the behavior of the TP-AGB models beyond R Hya’s classical constraints does not matter—with the production of a useful grid. The limit of 3 days’ run time is the main source of truncation for the fixed helium gird. The limit on max_number_retries is the main source of truncation for the helium-varied grid. Model summary files are generated at the conclusions of phases 1 and 2 as well as at the conclusion of phase 3 if the model is not terminated early. These are given names of the form TACHeB.mod (end of phase 1), AGB_seed.mod (end of phase 2), and AGB_terminal.mod (end of phase 3), respectively. Dividing the evolution into three parts allows us to relegate the asteroseismic calculations and use of velocity drag to phase 3 only. The AGB_seed.mod models, also made available with this grid, can also be used to start new calculations from the onset of the TP-AGB. GYRE is called once at each time step in phase 3.

As noted at the top of Section 3, we perform a sweep in two model input parameters only: initial mass and initial $Z$, with two different treatments of helium. We set $\alpha_{\text{MLT}}=1.931$ for all tracks. While we explored several values for the mass loss efficiency parameter $\eta$ initially, we did not add this as a dimension to the current grid (though it can be introduced later). The models we provide here cost just over 4.5 million CPU hours.

The calculations are performed as follows:

1. 1.

   We compare the TP-AGB (phase 3) component of each model to observational boundaries and isolate the regions where agreement is found. Agreement is strictly enforced, meaning a track must intersect the values and uncertainties prescribed in order to be considered in the next stages of analysis. Given the high variance among values quoted for the classical constraints in the literature, we use a relatively loose interpretation of the constraints listed in Table 1. The observational boundaries enforced in the analysis pipeline are summarized in Table 2.

   Deciding which models are candidates based on the theoretical periods requires making an assumption about the mode we are actually observing. We must therefore compare each model against the observational boundaries twice: once using the FM periods and once using O1 periods (for reasons discussed in Section 4, these are the only two possibilities.).

2. 2.

   Once the observationally consistent region of a TP-AGB pulse spectrum is identified, we use a signal processing routine from scipy.signal (Virtanen et al., 2020) to isolate the thermal pulses in this domain. Because thermal pulses contain sharp local maxima, this can be automated by specifying prominence and time step spacing criteria, but no one combination of these parameters is reliable for $\sim$3000 models, each of which contains anywhere from $\sim 5$ to $\geq 50$ pulses.

   Peaks that do not correspond to the period (luminosity) maximum within the power-down phase are sometimes identified, but most of these are discarded on the basis of either morphology or inconsistent physics: for example, if the hydrogen-to-helium burning ratio is not consistent with undergoing an active thermal pulse, set by $L_{\text{H}}/L_{\text{He}}<10$.

   In some cases, a secondary feature known as the helium sub-flash is identified, as demonstrated in Figure 7. As we have no physical reason to exclude sub-flash detections, they are are permitted. However, in the majority of cases the sub-flashes are excluded de facto due to their shallower slopes, enforced as follows.

   

   Because the more recent period observations have much denser sampling and much smaller uncertainties than the earliest measurements, these will inappropriately bias the pulse fit towards the last $\sim$third of the data unless counteracted manually. We use a “slope hardness” parameter, $s$, that sets a consistency threshold requiring

   $$\dot{P}\approx\Delta P/\Delta t\geq s/\Delta t$$

   (2)

   to avoid fitting pulses that are too shallow overall. We perform several tests and find that a setting of $s=75$ does the best job of identifying good fits and excluding poor ones. We caution, however, that the value chosen for $s$ can have a significant impact on the solution space. We computed our results for three settings of this parameter: $s=50$, $s=75$, and $s=100$. We focus on solution maps using $s=75$ henceforth, but results for all values are provided in Appendix B. Overall, this method identifies the correct regions of interest about 95% of the time.¹¹¹¹11Not perfectly— instances of incorrect pulse identification and indexing are documented in the data visualizer. These do not affect our results.

3. 3.

   For each pulse identified, the observational data are adjusted to center on the pulse, as shown in Figure 8. A data file that replaces the relative ages of the period observations with absolute ages from the models is generated and given a name beginning with auto_age_fit. Any given TP-AGB spectrum may host several acceptable pulse solutions (depending on how strictly observational agreement is enforced—see Appendix B) but typically there are between zero and $\sim 5$ auto_age_fits per model.

   As the age of R Hya is not well constrained (see Section 8.4), and realistic precisions on stellar age determinations are several orders of magnitude larger than the entire duration of a thermal pulse spectrum anyway (e.g. Tayar et al. 2022; Joyce et al. 2023), we treat absolute age as arbitrary but preserve the relative timing of the period measurements in our fits.

4. 4.

   Each period–age relation generated in step 3 is evaluated against its corresponding pulse on the MESA evolutionary track with the same mass and metallicity. In the immediate area of the pulse, the theoretical model is fit with four cubic splines: one curve relating each of period (FM or O1, depending on hypothesis), $T_{\text{eff}}$, luminosity, and radius to age.

   The cubic splines are sampled to generate theoretical $L$, $T_{\mathrm{eff}}$, $R$ and period values corresponding to each measurement. These are compared point-by-point to the observations using a weighted root-mean-square error (WRMSE) framework:

   $$\text{WRMSE}=\sqrt{\frac{\sum_{i=1}^{n}w_{i}\cdot(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{n}w_{i}}}$$

   (3)

   where:

   	$\displaystyle n\quad\text{is the number of observations or data points,}$
   	$\displaystyle y_{i}\quad\text{is the observed value for the }i\text{th data point,}$
   	$\displaystyle\hat{y}_{i}\quad\text{is the predicted value for the }i\text{th data point,}$
   	$\displaystyle w_{i}\quad\text{is the weight assigned to the }i\text{th data point.}$

   The weightings $w_{i}$ are given by

   $$w_{i}=\frac{1}{\sigma_{y,i}^{2}}$$

   where $\sigma_{y,i}$ is the uncertainty associated to observation $y_{i}$.

   When goodness-of-fit is determined according to the 210 d period measurements alone, we call the statistic $P_{\text{WRMSE}}$, or $P_{\text{W}}$.

5. 5.

   In the case where we wish to factor agreement with $L$, $T_{\mathrm{eff}}$, and $R$ into our consistency metric along with $P$, we compute Equation 3 for all quantities and combine the four WRMSE values two different ways. The first is a harmonic mean:

   $$H_{\text{WRMSE}}=\frac{4}{\frac{1}{x_{1}}+\frac{1}{x_{2}}+\frac{1}{x_{3}}+\frac{1}{x_{4}}}$$

   (4)

   where $x_{1},\ldots,x_{4}$ are the WRMSE values from (3) for $P$, $L$, $T_{\mathrm{eff}}$ and $R$. We abbreviate this as $H_{\text{W}}$.

   The second is a pseudo-$\chi^{2}$ score normalized such that the combination of all classical observational components is weighted equally with a seismic (period-only) score:

   $$S_{\text{WRMSE}}=\sqrt{\frac{1}{3}(L_{\text{W}}^{2}+T_{\text{W}}^{2}+R_{\text{W}}^{2})+P_{\text{W}}^{2}}$$

   (5)

   where $P_{\text{W}}$, $L_{\text{W}}$, $T_{\text{W}}$ and $R_{\text{W}}$ are the values from Equation 3, as above. We abbreviate this as $S_{\text{W}}$.

An interactive data visualization application is provided with this grid. We provide data files containing the calculations necessary to reproduce all of the realizations in Figures 9 and 10 and their counterparts for the helium-varied grid, as well as additional parameter variations explored in Appendix B.

Some of the visualizer’s capabilities include mouse-over information, pulse identification plots, and links to the evolutionary tracks associated to each model. These are demonstrated in Figure 24 in Appendix C. However, the documentation in the Github repository itself is most comprehensive: https://github.com/mjoyceGR/AGB_grid_visualizer.

Given that a substantial fraction of models were terminated based on run time limitations (72 or 120 hours maximum for the fixed helium or helium-varied grids, respectively) or convergence difficulties (max_number_retries $=500$) rather than by virtue of reaching their physically motivated stopping conditions, it is possible that any such model was stopped before encountering the thermal pulse that would provide the best fit to observations. Run incompleteness could therefore impact the score shown for a model in the heat map, and non-uniformity of incompleteness (e.g., models with very super-solar abundances are more likely to encounter convergence failure than solar-like models) could impact trends in the heat maps overall.

That said, premature model termination is not a concern for R Hya, which is typically best fit by one of the first 15 TPs regardless of helium variation. Where models do max out their run time, they have evolved far beyond the observational constraints for R Hya. This can be confirmed by inspecting the pulse identification figures available for every model in the grid visualizer.

There are, however, more significant limitations to this grid in terms of its utility as a general AGB fitting tool. Chief among these is the lack of variation in mass loss. The number of thermal pulses an AGB model undergoes is highly sensitive to the choice of mass loss efficiency, $\eta$, and mass loss prescription. The strength of the TPs and occurrence of third dredge-up depend on the mass of the envelope, and how the envelope mass changes over the course of the TP-AGB phase is likewise sensitive to these mass loss settings (see, e.g., Ventura & D’Antona 2005; Stancliffe & Jeffery 2007; Marigo et al. 2020; Karakas et al. 2022.)

Similarly, this grid does not consider variations in $\alpha_{\text{MLT}}$ or choice of opacity tables, the impacts of which have been documented in Cinquegrana & Joyce (2022b) and Cinquegrana et al. (2023). We do not explore variations in solar scale. We also note that the quantity $\Delta Y/\Delta Z$ remains fixed to 2.1 for the grid that scales $Y_{\text{init}}$ with $Z$. Not only is there considerable tension in the literature regarding the appropriate choice for $\Delta Y/\Delta Z$, but it may vary depending on the region of the Galaxy, local chemistry, number of enrichment sources (e.g., supernovae) in the vicinity, and time.

It is likely that scores would change in response to modifying any one of these stellar modeling assumptions, and possible (though less likely, based on explorations done here) that overall trends could shift as well. Such questions are worth exploring in future investigations. Nonetheless, this grid is an excellent resource for determining (preliminary) fundamental parameters for AGB stars and provides a much more accessible alternative to performing one’s own AGB calculations.

Finally, the grid calculations presented here are linear and adiabatic, and neither of these simplifications applies in earnest to AGB stars and their pulsation modes. We explore non-linear and non-adiabatic considerations further in Section 9.

The FM and O1 composite maps show very similar features in their solution spaces. Both maps feature a diagonal ridge tracing the preferred solutions, shown in red (warmer colors). The ridge in the FM map runs from around $1.9\,M_{\odot}$ up to almost $4.0\,M_{\odot}$ with increasing metallicity up to [Fe/H] = 0.7. Along this ridge, solutions in the upper region starting from $2.5\,M_{\odot}$ and nearly solar metallicity and upwards have lower scores than those below this. The preferred solutions in the O1 composite heat map (panel F of Fig. 10) trace a smaller area and form a more vertical ridge that trends towards lower masses. This ridge runs from $1.5\,M_{\odot}$ to $2.6\,M_{\odot}$, between [Fe/H] indices of $-0.6$ and $+0.8$. In both modes, the ridges are truncated where models cease to intersect the observational boundaries for any pulse.

The difference in slope between the ridges in Figures 9 and 10 reflects variation in the O1/FM period ratio, which depends on the physical parameters of the models. As we have shown in previous works on variable evolved stars, the strong dependence of the O1/FM period ratio on quantities like luminosity can be used to further constrain the model fits. This technique was used to characterize the TP-AGB star T UMi and the red supergiant Betelgeuse (Molnár et al., 2019; Joyce et al., 2020), but unfortunately cannot be applied to single-mode stars like R Hya.

Treating all solutions constituting the primary ridges in the composite maps (F panels, Figures 9 and 10) as equally valid, the best-fitting parameters of R Hya follow well-behaved quadratic scaling relations between mass and metallicity, depending on the mode assumed. For the FM, this relation is approximately

and for the O1,

These are shown overlaid on the heat maps in Figure 11.

Regardless of mode assumption, the best solutions fall towards the higher-mass edge of the regions permitted by the classical constraints. Based on the modeling results alone, we can already surmise that the FM is more likely than the O1: the number of good solutions is larger and covers a wider range of (sensible) initial masses and metallicities in this case. However, we can further discriminate among all good solutions by comparing against independent parameter estimates for R Hya in the literature and taking into account other sources of physical information, such as the presence of third dredge-up indicators, described in detail in Section 10.

Table 3 summarizes past attempts at determining the fundamental properties of R Hya, including mass, metallicity, and age. In cases where the parameter determination was made using AGB stellar models, the value is indicated in boldface. Figure 12 shows some of these overlaid on the FM and O1 composite ($S_{\text{W}}$) maps, with colors muted for visualization.

We first note that there is a huge range of masses and luminosities estimated for R Hya in the literature, as indicated also in Table 3, and very few estimates for its metallicity or age. The mass boundaries derived by Zijlstra et al. (2002) (similarly, Wood & Zarro 1981 and Eggen 1985) intersect with a subset of our best-fitting models between [Fe/H]${}\sim-1.0$ to $-0.5$ under the FM assumption, but considerably higher metallicities, [Fe/H]${}\sim 0$ to $+0.3$, under the O1 assumption. This is indicated by the blue hatched region in Figure 12.

The higher metallicities required by the O1 map are less likely, given R Hya’s location in the solar neighborhood, so the mass constraints suggested by most of the references in Table 3 provide further evidence that the period we have measured corresponds to the FM. One exception is the mass range reported by Decin et al. (2020) (green hatched region in Figure 12), which overlaps very little with our preferred FM solutions, barely intersecting them at the lowest masses and metallicities consistent with the observational boundaries.

The O1 solutions show slightly better agreement with Decin et al. (2020) around [Fe/H]${}=-0.5$, again for the lowest masses. Mass estimates from Fadeyev (2023) (red hatched region) are significantly higher than the observational boundaries permit in either case. Even assuming that the FM scaling relation (Equation 6) is valid above $4.0M_{\odot}$, a mass greater than $4.8M_{\odot}$ would require a metallicity exceeding the maximum metallicity of the Galaxy, indicated by the black shaded region at the top of both panels in Figure 12. Hence, our results are inconsistent with Fadeyev (2023).

Eggen (1985) suggested that R Hya might be a member of the Hyades moving group, which would constrain its age and mass to be about 0.5–1.0 Gyr and $\sim 2\,M_{\odot}$, respectively. This is the origin of their age value in Table 3. However, more recent results indicate that, while the Hyades cluster contains stars that also belong to the tidal streams of Hyades, the moving group itself originates from a resonance with the spiral arms of the Milky Way (Famaey et al., 2007; McMillan, 2011). According to Famaey et al. (2008), about half of the members are not coeval with Hyades, and hence this is not useful as an age constraint for R Hya.

While we cannot connect R Hya to any cluster directly, we can use various period–age relations to estimate its age. However, different age relations give very different answers, especially in the 300 to 500 d period range (Zhang & Sanders, 2023). Here, we assume that the 350 d period is a good approximation for the out-of-pulse period of R Hya, but calculate the ages for 300 and 500 d values as well, listing them in parentheses: In the 350 d case, the relation of Eggen (1998) gives an age of 6.9 Gyr (3.6–8.5 Gyr); Zhang & Sanders (2023) gives a very similar value of 6.5 Gyr (4.1–7.4 Gyr). In contrast, the relation of Grady et al. (2019) finds a much younger age of 1.01 Gyr (0.17–1.83 Gyr). The most detailed relations are presented by Trabucchi & Mowlavi (2022), who provide limits on the parameter space: in this case, the 350 d period corresponds to an age range of $0.5^{+2.6}_{-0.3}$ Gyr. Differences in the calibrating samples on which these period-age relations are based could be responsible for the large variance in the solutions they provide: any relation is inherently biased towards the cluster members from which it is calibrated and will therefore not be universally applicable to all Mira and TP-AGB stars.

Another complicating factor is that these calibrations assume a smooth evolution of pulsation period with age, but stars experience large, non-monotonic changes in period over the course of a TP. As stars evolve along the asymptotic giant branch, their radii and thus pulsation periods increase over time, on average. However, this overall trend of period increase with time is decorated with the high-amplitude variations of thermal pulses on much shorter timescales, meaning the star can enter the same period regime at very different points along its AGB evolution, leading to a period–age degeneracy (as discussed in Section 4 and shown in Figure 5). As such, period–age relations may not be the appropriate age inference method for stars suspected to be actively undergoing a TP.

However, our more sophisticated fitting method breaks the degeneracy present in simple period–age relations. Figure 13 shows the age of the best-fitting thermal pulse from each TP-AGB spectrum on the same mass–[Fe/H] plane as Figures 11 and 12. The fixed-helium assumption applies to the left-hand panels, and the varied-helium treatment applies to the right (with higher model failure rate, as discussed in Section 3, indicated by grey squares). Ages are in Myr, converted to $\log_{10}$ units and indicated by color. The thick red line in both panels is a broadened version of the quadratic fits shown in Figure 11 and indicates the ridge of best agreement for the FM (top) and O1 (bottom) solutions.

Our best-fitting models indicate a low age for R Hya: FM solutions have ages below 1.0 Gyr regardless of helium assumption. As shown in Figure 13, FM solutions trace a ridge ranging from 800 Myr (1 Gyr with helium varied) at the lowest masses to about 300 Myr at the highest, whereas overtone solutions give a somewhat higher age range of 1.0 to 1.6 Gyr (extending to 2.5 Gyr when helium is varied). These values are in good agreement with ages inferred from the period–age relations given by Grady et al. (2019) and Trabucchi & Mowlavi (2022), but not with the relations of Eggen (1998) and Zhang & Sanders (2023), which overestimate the ages significantly.

Trabucchi & Mowlavi (2022) also provide period–mass relations based on stellar models.¹⁴¹⁴14The lower mass boundaries are visual estimates based on Figure B.1 of Trabucchi & Mowlavi (2022), because the coefficients for the upper edge of the O-rich relation listed in their Table B.2 do not reproduce the line in the figure. Using the same set of periods from the previous section for comparison, we find initial masses of $2.72\,M_{\odot}$ (1.52–3.98 $M_{\odot}$), for the 350 d period. The mass range is slightly lower for a 300 d period (1.26–3.28 $M_{\odot}$) and considerably higher for 500 d (2.45–6.28 $M_{\odot}$).

An alternative way to estimate masses for TP-AGB stars is through the ¹⁷O/¹⁸O isotope ratio. According to the model calculations of De Nutte et al. (2017), this is a sensitive and reliable probe of stellar mass throughout the TP-AGB phase. Using the isotope ratios measured by Hinkle et al. (2016), Decin et al. (2020) found a significantly lower mass of $\sim 1.35\,M_{\odot}$. Taking the uncertainties in the isotope ratio into account, a mass range of 1.0–1.5 $M_{\odot}$ can be estimated. We note, however, that Hinkle et al. (2016) found a ratio of $1.0\pm 0.4$ for $\chi$ Cyg, whereas De Nutte et al. (2017) cites a ratio of $2.0\pm 0.5$ for the same star, so the uncertainty in this method or in the isotope ratio could be even higher.

We list these estimates in Table 3 and mark them on Fig. 12. As the maps show, Decin et al. (2020)’s estimate based on the ¹⁷O/¹⁸O ratio clearly underestimates the mass of the star compared to the other solutions presented here: while the O1 ridge extends to masses that partially intersect this range, the favored mass ranges are significantly higher for both mode assumptions.

On the other end of the mass scale, the mass estimated by Fadeyev (2023) is strongly disfavored by our models. This can be explained by the fact that Fadeyev (2023) assumed that the earlier “hovering” (standstill) in period corresponded to the end of the power-down phase, whereas recent observations show otherwise.

The generic assumption of a 2 $M_{\odot}$ model made by Wood & Zarro (1981), Eggen (1985) and Zijlstra et al. (2002) falls comfortably within the FM and O1 ridges, but only at relatively low, sub-solar metallicities ([Fe/H] $\leq-0.4$) for the FM models.

The masses we calculated from the period–mass relations of Trabucchi & Mowlavi (2022) constitute a broad range that encompasses our best-fitting masses. We find that assuming 350 d for the fundamental mode gives the most consistent mass values from this relation: 1.52–3.98 $M_{\odot}$ (best fit at 2.72 $M_{\odot}$), suggesting that the out-of-pulse “equilibrium” period for the star is close to 350 d. The picture is more complicated for the first-overtone masses, since the mass–period and mass–age relations are calibrated by assuming fundamental-mode pulsation.

Homan et al. (2021) assume that R Hya is a multiple system, and the mass they provide is an estimate of the present-day mass of a central reservoir containing two components. While their reported mass range is comfortably consistent with many of our preferred models under either mode assumption, the appropriate comparison for Homan et al. (2021)’s results would be to the end state of binary evolution calculations, which we do not compute.

As shown in Section 2.5, the metallicity of R Hya is possibly super-solar. Metallicities in the disk of the Milky Way extend to about [Fe/H] = 0.5–0.52, based on results from various large-scale spectroscopic surveys such as APOGEE (Hayden et al., 2015), GALAH (Sahlholdt et al., 2022) and GES (Gaia-ESO Survey, Randich et al., 2022). Even higher metallicities have been claimed in the Galactic bulge (Hill et al., 2011; Bensby et al., 2013), but the kinematics suggest that R Hya is a disk population star. To avoid folding assumptions about R Hya’s membership of particular populations into our analysis, we treat R Hya’s metallicity agnostically and explore a wide range of values centered on the solar [Fe/H], ranging from $-1.0$ to $+1.0$. The latter is the upper limit stars can theoretically reach, according to the calculations of Cheng & Loeb (2023).

For the FM assumption, we find good solutions at almost every value of [Fe/H] ($Z_{\text{init}}$) shown in Figures 9 through 12, with the exception of the two lowest metallicities shown. However, calculations were also performed for [Fe/H]${}=-1.53$ and [Fe/H]${}=-2.2$ (not shown) in an exploratory grid (see Section 3), and no model with either of these metallicities, regardless of mass, intersected any of the observational constraints. They were thus excluded, and no further efforts to study metallicities between [Fe/H]${}=-2.2$ and [Fe/H]${}=-1.2$ were made on the basis of lack of observational consistency.

An “exclusion zone” set by an upper limit of $+0.52$ dex is indicated in Fig. 12. In the static helium case (left column of Figure 12), the FM solutions extend into this regime for masses above $\sim 4.0M_{\odot}$, suggesting that such masses are not realistic. Perhaps unsurprisingly, such models generally fail due to convergence difficulties in the helium-varied grid (right panels of Figure 12) when computing compositions with metallicities above this limit, as helium becomes nearly 50% of the star’s total composition.

The O1 metallicity range is narrower, beginning at $-0.5$ dex for the lowest masses in the fixed-helium grid but extending well beyond the Galactic upper limit, encompassing models with metallicities as high as $+0.9$ dex. The upper third of the O1 solution space is thus excluded on this basis, ruling out masses $\geq 2.6M_{\odot}$.

In addition to being more physically correct representations of the conditions of AGB stars, non-adiabatic models can calculate the linear growth rates of modes, indicating which modes are unstable to small perturbations. Mode stability for red giant and AGB stars has been investigated by many authors (see, e.g., Xiong & Deng, 2007; Takayama et al., 2013; Xiong et al., 2018; Trabucchi et al., 2019). As we touched upon in Section 4, in principle, non-adiabatic models could more definitively constrain R Hya’s pulsation mode. However, mode selection in the non-linear regime (i.e., observing which modes grow to full amplitude) and mode identification based on observations pose theoretical and practical difficulties. As Paper I demonstrated for T UMi, stars might go through mode transitions over the course of a TP, including entering a double-mode state. Mapping the changes in the stability parameter over the entire grid would increase the computational requirements significantly. Theref ore we elected not to constrain our grids with a growth rate criterion, but rather to map the full parameter space both for the FM and O1 modes.

Given that non-adiabatic solutions may differ from adiabatic ones for stars like R Hya, we recomputed a handful of our preferred models (described in more detail in Section 11) using GYRE in non-adiabatic mode. The results of one representative model ($M=2.9,Z=0.0216,Y=0.294$) are shown in Figure 14. We note here that the difference in periods predicted from running GYRE in adiabatic (magenta/pink) mode vs non-adiabatic (black/grey) mode is not significant for the (FM) pulses that are observationally consistent with R Hya. The discrepancy widens as the star evolves, becoming more significant at later pulses, as the star’s outer envelope becomes more stratified and the conditions become more strongly non-adiabatic. However, the difference does not, for any of the representative models explored, come close to exceeding the 60-d uncertainty already incorporated into the numerical experiment. The adiabatic vs non-adiabatic discrepancy is roughly t he same for O1 as for FM periods.

Based on this marginal discrepancy of, typically, $\sim$20 days or fewer, we conclude that the conditions in models that fit R Hya are only weakly non-adiabatic, and that the adiabatic approximation is valid for this particular star. Though the validity of this approximation does not extend to every pulse in every model over the entire grid, the substantial increase in computation time and structural resolution required to run GYRE in non-adiabatic mode is prohibitive. We therefore do not attempt to run a comprehensive grid in non-adiabatic mode, but we do provide the GYRE control file (gyre_non-ad.in) used to perform non-adiabatic calculations in the repository. We leave a more thorough description of these settings to the documentation on GitHub, but note that our non-adiabatic search uses the adiabatic eigenfrequencies as the trial roots for the non-adiabatic case (as documented in description of the ‘AD’ method for non-adiabatic calculations in the GYRE documentation). In this way, the adiabatic solutions presented in our grid still have utility regardless of how well the adiabatic approximation applies.

A limited alternative to running non-adiabatic calculations directly is provided by Xiong & Deng (2007) for the RGB. However, we found that their analytic approximation was (strongly) invalid for the AGB. Xiong & Deng (2007)’s stated domain of validity for the pulsation constant, $\log_{10}Q$, in their approximation is $\log_{10}Q=-1.30$ to $-0.9$, whereas calculated values of $\log_{10}Q$ for our models ranged from $0$ to $0.25$ along the AGB.

As discussed in Paper I, working with linear periods introduces uncertainties in our model fits. While small-amplitude variations in less luminous stars can be approximated well with linear periods, large-amplitude pulsations rearrange the envelope structure of the star, shifting the observed, non-linear periods away from the linear values. However, the degree, or indeed even the direction, of this period shift has not been conclusively determined, and results vary depending on the non-linear pulsation models different authors use. For example, Ya’Ari & Tuchman (1996) and Lebzelter & Wood (2005) found that the periods shorten in the non-linear regime, whereas periods lengthened in the models of Kamath et al. (2010) and Ireland et al. (2011). Dependence on initial mass and convective model prescriptions have also been documented Olivier & Wood (2005). While MESA itself has a non-linear pulsation module included, the RSP (Radial Stellar Pulsations; Smolec 2016) code was developed to reproduce less luminous RR Lyrae and Cepheid-type variations. MESA-RSP runs into dynamical instabilities at around 1000 solar luminosities, which is an order magnitude below the luminosity of R Hya.

Besides the uncertainties in various parameter choices, the main drawback of non-linear models is that they are very computationally expensive. Because of this, non-linear studies typically rely on only a small set of models, making extrapolations to other physical parameter regimes difficult. More recently, however, Trabucchi et al. (2021) calculated a denser grid of models that made it possible to deduce scaling relations for fundamental mode period differences between linear and non-linear models. For these fundamental-mode pulsations, Trabucchi et al. (2021) found that non-linear periods begin to shorten relative to their linear counterparts above a break point in radius (or in period), and eventually saturate at a constant maximum period for the largest stars. This phenomenon would undoubtedly affect our models, in effect lowering peak periods and the hardness of the power down phase, and longer periods will be affected more. For overtone solutions, Trabucchi et al. (2021) found these differences to be negligible, therefore our conclusions for the first overtone model grid are not affected by non-linear shifts.

To investigate the severity of the impact on our results, we implemented the scaling relations of Trabucchi et al. (2021) into our post-processing algorithms and re-ran our analysis using the rescaled period evolution data in lieu of periods computed directly with GYRE. However, since this scaling also introduces many implicit assumptions that are specific to the particular non-linear models used by Trabucchi et al. (2021) (and there is tension in the literature regarding the direction and magnitude of linear vs non-linear trends), we consider these scaled results only as a comparison to the linear model grid.

In Paper I and in the present study, we treated period constraints loosely, both to accommodate possible non-linear shifts and to account for the occurrence of TPs at discrete period values in the models that may or may not match the observed periods perfectly. However, in the example provided in Figure 15, we see that the discrepancy between linear and non-linear periods is roughly 70 days, which is slightly greater than our uncertainty allowance. The non-linear periods therefore differ enough from the linear calculations that it is worth investigating the impact on the global solution space. A side-by-side, linear vs non-linear comparison for the fundamental mode solution space is shown in Figure 16.

Globally, our results for R Hya do not change significantly when using the non-linear scaling relations in lieu of direct calculations with GYRE. It is worth noting that using scaling relations, which predict non-linear periods as a function of present-day mass and radius only, is less computationally expensive than running GYRE. However, we reiterate that these scaling relations are based on case-specific calibrations from one particular non-linear code, and the assumptions adopted in Trabucchi et al. (2021) may or may not be valid across our entire parameter regime. Rather than repeating the analysis of Section 8 for the right-hand panel of Figure 16, it is sufficient to observe that (1) using the non-linear scaling relations does have the anticipated impact of pushing R Hya’s best-fitting models towards lower peak periods—especially for lower masses and metallicities; and (2) it becomes more worthwhile to compute or extrapolate from non-linear models the more evolved (i.e. later in the pulse spectrum) the model is.

We proceed with our analysis using the linear and adiabatic calculations as the de facto results, but explore the impact of non-linearity on best-fitting pulse indices in Section 11.