© The Authors 2023

1 Introduction

Mira variables are located far into the low-temperature side of the Hertzsprung-Russell (H-R) diagram. In this region, Mira (Omicron Ceti) was the first star to have its pulsation period identified and, together with other Mira variables, was amongst the first periodic red variable stars to have been determined (see Hoffleit 1997, for historical remarks). Mira-type stars are extremely evolved stars at the tip of the asymptotic giant branch (AGB). Low-to-intermediate mass stars (with initial masses in the range 0.8 M_⊙ < M_* < 8 M_⊙) encounters the AGB phase during their stellar evolution. Variability is common among cool giant stars. Studies (see, for example, Kiss et al. 1999; Cioni et al. 2001; Soszyński et al. 2009) have revealed that the differences, in terms of luminosity amplitude and period regularity, are great.

The pulsation periods of AGB stars range from 80 to 1000 days. The regularity of the variation of the stellar luminosity over time allows the subdivision of the pulsating variable types within the AGB: Miras, semi-regular (SR), semi-regular variables with regular periodicity (SRa), and semi-regular variables with irregular periodicity, which can be subdivided depending on the presence (SRb), or lack of (SRc), any transient periods being observed (Becker 1998). Red supergiants (RSGs) are often classified as SRb and SRc stars, showing visual variability with an ambiguous period or amplitude. The pulsating AGB and RSG stars are collectively called long period variables (LPVs); however, Miras are distinguished by high amplitude luminosity variations (∆m > 2.5 mag), with a long and singular dominant period, in the low-order radial pulsation mode. Stars belonging to the SR, SRa, SRb, and SRc classes show smaller amplitude variations with multiple pulsation periods. Consequently, these properties have allowed the AGB stars to be identified through observation, in particular by the period-luminosity (P-L) relations. Miras follow a well-constrained linear P-L relation (Whitelock et al. 2008, 2009) corresponding to the fundamental mode branch, whereas the semi-regular variables have P-L sequences parallel to that of Miras and the positions of the sequences are thought to be dependent on which pulsation mode is dominant (Wood et al. 1999; Soszyński et al. 2013; Wood 2015).

Pulsations are not limited to stars that are highly evolved, but are common throughout the H-R diagram. Solar-like pulsations in the range of micro-magnitudes are due to stochastically excited oscillations, whereas Cepheids and RR Lyraes have higher pulsation amplitudes generated by the κ mechanism (Joshi & Joshi 2015). However, AGB stars differentiate themselves from such well-known types of pulsating stars in two important respects. First, their visual amplitudes are high. Second, a degree of variability is inherent in the pulsations. Considering these two properties together with their stellar parameters suggests the excitation mechanism of the pulsations is different than stochastic excitation or the κ mechanism (see relevant review chapter in Catelan & Smith 2015). Modelling considerations for the excitation mechanism are complicated, since many processes interact with each other. The cool red giants have extremely deep convection zones with large convective cells, often introducing surface inhomogeneities. The large convective cells actively being cycled in the interior, through a series of up- and downdrafts, mean there is a complex interplay between convection and pulsations, therefore non-linear models are needed for the high-amplitude pulsations.

Using 1D linear pulsation models, Trabucchi et al. (2019) were able to reproduce observed periods of overtone modes in evolved giants, but their predictions for the fundamental mode periods do not agree with observed relations. In recent work with non-linear models, however, Trabucchi et al. (2021) obtained fundamental mode periods more closely in agreement with observations. This indicates that convection and pulsations are not accurately described in a linear fashion; the interplay needs to be treated non-linearly. Still, the produced models are 1D in nature and the description of convection is parameterised according to the mixing-length theory; which, despite its usefulness, is insufficient to reconstruct accurate stellar interior temperatures and pressure structures (see reviews by Houdek & Dupret 2015; Kupka & Muthsam 2017, for a thorough discussion). The fine-tuning of the mixing-length parameter to reproduce such structures compliant with the laws of physics risks oversimplifying the physical processes active in stellar interiors. The need to use 3D modelling is exemplified by the fact that convection, turbulence and the excited oscillations are intrinsically 3D in nature. Utilisation of 1D models can complement 3D models in furthering our understanding of pulsations in evolved red giants; they can explore a wider range of stellar parameters and dynamics, being appreciably less computationally expensive than complete 3D models. A further summary on the use of 1D models of convection and pulsations is available in the introduction section of Freytag et al. (2017). In conclusion, it is hard to overstate the need for global 3D models in modelling pulsations of evolved red giants. CO5BOLD is a radiation-hydrodynamics code that numerically simulates global stellar convection. In simulations of cool red giants, self-excited pulsations exist as standing pressure waves. The developed surface convection cells are consistent or comparable with surface inhomogeneities seen from observations (Freytag & Höfner 2008; Freytag et al. 2017; Kravchenko et al. 2018; Paladini et al. 2018; Chiavassa et al. 2020).

In this work, pulsation properties of global 3D models of AGB and RSG stars produced with the CO5BOLD code were extracted and correlations against stellar parameters were investigated. In essence, the results of the simulations were treated as observed data to identify any significant physical trends. The paper is presented in the following manner. The setup of the 3D models, the model properties and how the pulsation properties were extracted are presented in Sect. 2. Section 3 outlines our main results, with relevant discussions presented in Sect. 4. Section 5 contains our final remarks.

2 Methods

In this section, we first give a short overview of the basic properties of the CO5BOLD models that we analyse in this paper. Then we discuss the methods used for extracting pulsation properties. In Sects. 2.2 and 2.3, the model st26gm07n001 is referenced as the representative standard model, since this model produced reliable pulsations qualitatively similar with most of our models. This model is a M_* = 1 M_⊙ mass AGB model with solar abundances, having the time averaged stellar parameters radius, luminosity, and effective temperature, respectively: R_* = 360 R_⊙, L_* = 7000 L_⊙, and T_eff = 2700 K.

2.1 Setup of 3D models

For our simulations, we used the CO5BOLD code (Freytag et al. 2012; Freytag 2013, 2017), that numerically integrates the coupled non-linear equations of hydrodynamics and radiation transport. Compressible hydrodynamics equations allow the modelling of travelling pressure waves, standing acoustic modes, and transsonic convective flows in the stellar interior, as well as shocks in the atmosphere. A short-characteristics scheme solves the equations of non-local radiative energy transport, appropriate for the optically thin atmosphere and the optically thick interior. The gas opacities are tabulated and are the result of a merger between tables generated from stellar-atmosphere codes such as PHOENIX (Hauschildt et al. 1997) or COMA (see Aringer 2000; Aringer et al. 2016) for the outer layers and OPAL data for the stellar interior (Iglesias et al. 1992). No radiation pressure on gas is included, so far.

The tabulated equation of state accounts for the ionisation of H and He, the formation of H₂ molecules, and a representative neutral metal. Tables that take the frequency dependence of opacities into account are necessary to model the temperature structure of the outer atmospheric layers correctly. That is crucial, for example, for the formation of dust grains (Höfner & Freytag 2019). However, most of the models analysed here only use grey (Rosseland mean) opacity tables. These are sufficient for the study of stellar pulsations, and save significant computer time. As the radiation-transport with grey opacities already dominates the overall run time of a simulation, it slows down with frequency-dependent tables by a factor of about two (because of the further reduced radiative time step) times the number of opacity bins employed – typically five or three for global models (see Chiavassa et al. 2011; Höfner & Freytag 2019, respectively).

All models use a 3D cubical computational box and a Cartesian grid, with a grid dimension ranging from 171³ to 765³ cells (see Table A.1). For most models the computational box comprises the star and has equidistant grid axes in all directions. For a few models, this domain was augmented with an outer box, where the distance between the grid points gradually increases with the distance from the stellar centre. That allows the tracking of dust clouds (see Höfner & Freytag 2019, 2022) but is irrelevant for the properties of the pulsating star. The gravitational potential is a prescribed formula assuming that most of the mass is located in a compact, not resolvable stellar core. It follows an 1/r law in the outer layers and is smoothed in a core region with radius r₀ (typically 20–25 % of the stellar radius). While several formulas have been tried, the most often used is of the form $$\varphi \left(r\right)=-G*M*{\left(r^4+r_0^4\right)}^{-1/4}.$$(1)

Inside a sphere with radius r₀, energy is injected at a constant rate according to the desired luminosity. A drag force is applied that reduces the local velocity in each grid cell in the core region by a small amount per time step, independently of the size of the flow structures. This is in contrast to the action of an artificial viscosity, that is sensitive to velocity gradients. There are options, for instance, to make the drag force stronger in the very centre than in the outer layers of the core. The drag force, acting only in the small core region, damps the dipole flow globally. It also damps smaller-scale convective flows in this core region, where the gravitational potential is smoothed, anyway, and we cannot expect a good representation of the conditions by our simulations. Due to the small amplitude of radial pulsations in the core region, we expect only a small effect of the damping on this type of oscillations. In effect, the drag force leads to a smooth decline of the velocities from the boundary of the core region towards the very centre, with a rate depending on the detailed settings.

The initial models for the first global simulations were produced by integrating 1D hydrostatic stratifications and distributing the data onto a 3D sphere. This led to long transition phases, in which convective structures developed first in the surface layers (see Fig. 1 in Chiavassa et al. 2011) and only later extended further down. Such a procedure works well when initial 1D and final 3D stratification are similar, for example, when the stratification is close to radiative equilibrium or nearly adiabatic and the dynamical pressure due to convection and pulsations is negligible. However, in other cases it is preferred to start from an existing well relaxed model. Straight-forward transformations, for instance, the interpolation to a different grid or the expansion to a larger computational box are done via external scripts. In addition, some control parameters such as stellar mass or core radius can be changed gradually over the first few pulsation cycles during a simulation. During this phase, the density and the internal energy can be adjusted in all grid cells in each time step. That allows the smooth transition from one set of stellar parameters to another, shortening a long transition phase except for the scaling time interval itself. We monitor, for instance, the emerging luminosity, the total energy, entropy, pressure, and temperature in the core region, and the radius of the star over time to ensure our models are relaxed. The initial transition phase is omitted from further analyses.

Our dataset shown in Table A.1 comprises models with current stellar masses of 1 and 1.5 solar masses, representing AGB stars, as well as models between 5 and 12 solar masses, that we take as representative for RSG stars. Earliest CO5BOLD models are omitted, because often the time sequences are rather short for a reliable determination of the pulsation properties, and the simulations were replaced by newer ones with higher numerical resolution spanning longer time intervals. However, we note that the models analysed here are computed with different versions of CO5BOLD, with unequal grid sizes, and other differences in the numerical setup to maximise the available range of stellar parameters.

2.2 Definition of the stellar radius

In contrast to other stellar parameters, which are more or less explicitly set in the simulations, the stellar radius, R_*, is a resulting property of the models. It is important to have a reliable and consistent method for determining the stellar radius of the model. Not only would the pulsation extraction process be affected by an inexact radius, but it would also risk misrepresenting the stellar parameters. To define the stellar radius, one method we used was the same as Freytag et al. (2017), who relied on the average of the luminosity and temperature profile of the star, and solved for the radius following the Stefan-Boltzmann law 4πσR²T(R)⁴ = L. Although this definition works well in most cases, it can be unreliable when spherical averages are applied in the presence of significant opacity fluctuations arising due to variations in the surface density and temperature. For example, in models with similar stellar parameters that only differ in whether they are based on grey or non-grey radiative transfer, the radii derived from the latter are calculated to be greater. This is primarily due to fluctuations in the density that naturally occur in such simulations. The fluctuations may introduce sudden opaque structures that can be misinterpreted as the outer boundary of the photosphere which cause the gradient of the temperature profile in the surrounding atmosphere to be less well-defined. Hence definitions that rely on the explicit detection of temperature changes in the atmosphere, including averaging of the temperatures, can lead to inaccurate values of the radius.

Six further values for the radius have been used in this study. Four were defined by setting the requirement of the Rosseland mean opacity (τ_Ross) to equal 2/3, 1, or 10, or at the first instance of a spherically averaged temperature of T_avg = 8000 K going outwards from the centre of the star. Whilst using τ_Ross = 2/3 or 1 has practical and physical applications in relation to the photospheric radius of a star, the choice of τ_Ross = 10 was arbitrary in order to probe regions where it is substantially optically thick, therefore representing an estimate of the inner radius of the star. A similar motivation was behind defining a radius at T_avg = 8000 K. The temperature of 8000 K was chosen since it consistently captures the region where the inner temperature profile experienced the steepest decline going outwards from the centre of the star. This distinct region signified a transition region which accordingly may be taken as the radius of the star.

A new definition is devised in the work presented here, where the radius was determined by the point of the first local minimum of the entropy after a spherically averaged temperature of T_avg = 8000 K from the centre was reached. The entropy minimum marks the thin transition region between a convectively unstable stellar interior and a convectively stable stellar photosphere, and characterises the radius of the star in linear stability analyses and mixing-length theory (see Abbett et al. 1997; Trampedach et al. 2014; Magic et al. 2015). For the standard model, the spherically averaged entropy variation across the star is shown in the left panel of Fig. 1, where a local minimum can be clearly seen. The location of the entropy minimum depends on the pulsation phase. On the right panel of Fig. 1, the radial coordinate of the entropy minimum was plotted as a histogram and illustrates the contraction and expansion phases of the radius, therefore the minimum and maximum radial coordinates respectively, due to pulsations.

We set an additional definition for the radius for the purpose of capturing the transition region between standing and travelling waves. As will be shown in Sect. 3.5.1, this transition region is helpful in verifying which option for the stellar radius definition is best. Capturing the transition region involved determining the extent to which the pulsations are in phase throughout the layers near (just below, in, and just above) the stellar radius. This determination can be broken down into five steps: (1) the phase of the averaged radial velocity, ϕ_mean, over a thin height across the radius (following the entropy minimum definition) was computed; (2) the phase of the radial velocity of the layers 30% below and above the radius was computed and the phase difference, Δϕ, relative to ϕ_mean was calculated; (3) the standard deviation of Δϕ across the region was multiplied by 2 and used as a threshold; (4) new lower and upper boundaries of the radius were determined from the intersection between Δϕ and the threshold; (5) averaging was done between the lower and upper boundaries to represent the new radius, determined by the phase.

In total, we consider seven definitions for the stellar radius; each definition can provide different insights into the models and has its own advantages and disadvantages. As described in detail in Sect. 3.5, using the entropy minimum definition provides the best radius definition in the context of analysing the pulsation properties.

Fig. 1 Entropy profile of the standard model st26gm07n001. Left panel: temporal mean value of the entropy as a function of the radial coordinates, with the shaded area being the corresponding standard deviation and the zoomed-in axes shows the first local entropy minimum after the first instance of a spherically averaged temperature of 8000 K away from the centre of the star. Right panel: count values of where the entropy minimum was located over a simulation (in bins of size 2 R⊙).

Fig. 2 Diagnostics of the spherically averaged radial velocity for the standard model st26gm07n001. Top left: radial velocity as a function of time and radial coordinates, where the solid and dashed lines mark the radius and 80% of the radius, respectively. Bottom left: average and standard deviation of the radial velocity (in blue and shaded in grey, respectively) within the layer of 80–100% of the radius. Top right: FFT power spectrum (values are clipped to the maximum power calculated in between the 80–100% of the radius) as a function of the frequency and radial coordinates. Bottom right: averaged and normalised FFT power spectrum across the aforementioned layer. Here the shaded region is the extent over which the weighted averaging (see text) to get a characteristic pulsation frequency was done, while the solid and dashed black vertical lines indicate the computed characteristic frequency and its spread, respectively.

2.3 Pulsation period extraction process

As outlined in Freytag et al. (2017), to solely analyse radial motions, the radial mass flux and mass density were averaged over spherical shells (spherically averaged quantities are denoted with <.>_Ω, with Ω being the solid angle of a sphere). The ratio of the averaged radial mass flux, <ρv_radial>_Ω(r, t) and mass density, <ρ>_Ω(r, t), is taken as the radial velocity, thus becoming a function of radial distance and time <v_radial>(r, t). The corresponding spherically averaged radial velocity for our standard model as functions of time and radial coordinate are as shown in the relevant panels of Figs. 2 and 3. There are three regions of interest in the top left panel of Fig. 2: (1) the deep inner core shows irregular radial velocity structures, which arise from fluctuations in the convective flow; (2) around the photospheric layer, periodic variations of inward and outward flows of the radial velocity suggest regular periodic movements of matter due to pulsations; (3) above and beyond the outer atmosphere, where travelling shock waves dominate. The spherically averaged radial velocity under-represents the actual radial velocity amplitudes. Observed motions, derived from high-resolution spectroscopy line profiles, around the photosphere of Mira variables reach about 10 km s⁻¹ outwards, and about 15 km s⁻¹ inwards (see Hinkle 1978; Hinkle et al. 1984). The radial velocity calculated at each coordinate point for a single snapshot in the computational box, as presented in the right panel of Fig. 3, shows close agreement to the amplitude of the observed motions.

The fast Fourier transform (FFT) was applied to the radial velocity, to transform the signal into the frequency domain. Before the FFT, three pre-processing steps were applied to the input signal, which led to the most consistent results for the resulting FFT: (1) subtraction of the temporal average of the signal from the signal itself; (2) tapering with the Tukey window, which decreases the amplitude gradually at the start and end of the signal (thereby helps with dealing with possible discontinuities between the beginning and end of the signal); and (3) zero-padding to increase the length of the input signal by a factor of 10, (adding zeros at end of the signal produces higher resolution as there are more frequency points in the FFT). All the models presented in this work have sufficiently small sampling interval in the temporal domain to ensure a high FFT sampling rate thus fulfilling the Nyquist criterion safely.

The power spectrum at each spatial coordinate is shown in Fig. 2 and the region with the highest intensity can be taken as the first order estimation of the pulsation frequency. To improve on this estimation, a thin region below was selected, and the power spectrum was averaged over the layers in this region. The chosen region had to be selected strategically: too high above R_*, the shocks would dominate the power spectrum and the frequency extracted would not be the pulsation frequency; on the other hand, too low below R_*, convective flows may dominate the power spectrum since the pulsation amplitude decreases significantly, therefore obscuring the dominant pulsation frequency. With these effects in mind and after a calibration study, the optimum region to perform the FFT analysis was calculated to be in a layer between 80–100% of the radius of the star. The final height-averaged and normalised power spectrum is shown in its corresponding panel in Fig. 2, where dominant peaks are observable.

The ideal power spectrum would be primarily dominated by a single peak but, due to fluctuations to the pulsation amplitude and possible phase changes, multiple peaks may surround a dominant peak. The peaks within a region surrounding the main peak should be considered as part of the pulsation mode. To account for the smaller peaks, lower and upper bounds of the region to consider were set to 40% below and above the frequency of the strongest peak, respectively. Weighted averaging within the width was done to derive the mean frequency of the power spectrum. The weighted mean of the frequency within the interval was computed as $$f_{\text{mean}}=\frac{{\displaystyle {\sum }_{\text{i}=0}^{I}P\left(f_{\text{i}}\right)f_{\text{i}}}}{{\displaystyle {\sum }_{\text{i}=0}^{I}P\left(f_{\text{i}}\right)}},$$(2)

where I is the number of all the frequency points within the width (which varies with the time resolution of each model) and P is the corresponding power of the FFT power spectrum at the specific frequency point f_i. The corresponding weighted root-mean-square (RMS) is then $$f_{\text{RMS}}={\left[\frac{{\displaystyle {\sum }_{\text{i}=0}^{I}P\left(f_{\text{i}}\right){\left(f_{\text{i}}-f_{\text{mean}}\right)}^2}}{{\displaystyle {\sum }_{\text{i}=0}^{I}P\left(f_{\text{i}}\right)}}\right]}^{1/2}.$$(3)

The characteristic pulsation frequency and spread were determined from the values f_mean and f_mean ± f_RMS, respectively, and fit into the FFT power spectrum as shown in the vertical lines in the bottom right panel of Fig. 2. Due to occasions where a considerable number of smaller peaks were present outside the relative width, an adaptive width was utilised. The adaptive width extends the initial width only if the outer peaks contributed to a change in the final f_mean by at least 2%. This adaptive width did not extend indefinitely as it was also restricted within a range of extension. This is to avoid accounting for higher harmonics due the pulsations or possible overtone modes (where both would occur around twice the initial f_mean). Thus the adaptive width determination was useful to ensure all or most of the power due to the pulsations was contained within the width used in the weighted averaging. Finally, to convert into the time-domain, the derived pulsation period is simply P_puls = 1/f_mean, and its spread σ_puls = P_puls (f_RMS/f_mean). Since there were no nodes between the centre and surface of the star in the radial velocity profiles (evidenced in Figs. 2 and 3) and no appreciable power was present below the frequency of f_mean in the FFT space, we take P_puls to be the pulsation period of the fundamental mode.

Fig. 3 Radial velocity profiles for the standard model st26gm07n001. Left: snapshots of the spherically averaged radial velocity as a function of the radial distance, covering 4 yr of simulated time (about three pulsation cycles). Right: radial velocity computed at each cell of the computational box determined from a snapshot at time t = 14.2 [yr], presented to show the number of points (N) having the radial velocity at a certain distance from the centre of the box. The vertical lines mark the location of the stellar radius.

Fig. 4 Comparison between the different width methods, relative (Rel) or adaptive (Ada) in superscript, in the derived pulsation periods and spreads (left and right panels, respectively). Two models that have extremely extended atmospheres were omitted from the right panel to preserve clarity of the plot (as the values are around 0.8 [yr]).

3 Results

Given the self-excited nature of the pulsations in the models arising from pressure sound waves due to convection throughout the star (see Freytag et al. 2017, for additional background), the pulsation properties and correlations against various stellar parameters were investigated in this work.

3.1 Extracted pulsation periods and their spreads

The pulsation periods derived using either the relative or adaptive width method are compared in Fig. 4. Clearly, there is little scatter in the derived periods. There is a bigger scatter in the spread, where the adaptive method more consistently results in higher spreads in the pulsation period. This can be expected when the non-linearity of the interplay between pulsations and convection becomes considerable, as there are more peaks surrounding the dominant signal. With more appreciable peaks and the broader width being used in the weighted averaging procedure, larger values on the spread of the pulsation periods are expected when using the adaptive width; which is also suggestive of strong interaction between the pulsations and convection. The overall trend, therefore, of higher divergence from the one-to-one ratio at higher values in the σ_puls plot is reasonable. While it is important to keep the ${\sigma }_{\text{puls}}^{\text{Ada}}-{\sigma }_{\text{puls}}^{\text{Rel}}$ plot in mind, given the good correlation between the two methods of calculating the P_puls, unless otherwise stated, the results presented from now on are using the pulsation properties derived with the relative width method.

3.2 Representative examples

This subsection details the information presented in Figs. 5 and 6, where each column corresponds to a model with a different mass (1.5 M_⊙, 1 M_⊙, and 8 M_⊙). In Fig. 5, the bolometric luminosity variation and the averaged radial velocity in the layer between 80–100% of R_* are presented in the first and third rows, respectively, to indicate the amplitude of the pulsations in the models. The bolometric luminosity profiles of the three models show semi-regular variability (of the type SRa). The amplitudes of the 1.5 M_⊙ and 1 M_⊙ AGB models are low compared to those of typical Miras (see Feast et al. 1982; Whitelock et al. 2000). Disregarding any extinction (circumstellar or interstellar), the two AGB models show apparent bolometric amplitudes of around ∆m_bol = 0.14 mag, in contrast to most of the Mira-like pulsations observed in Whitelock et al. (2000), which had inferred amplitudes of higher than ∆m_bol = 0.5 mag (after bolometric corrections were applied).

The spherically averaged radial velocity and the corresponding averaged power spectrum after performing the FFT procedure outlined in Sect. 2.3 are shown in the second row of Fig. 5 and Fig. 6, respectively. The fourth row in Fig. 5 shows the recovered radial velocity signal in the time-domain, from an inverse FFT after masking was done in the frequency space to isolate only the pulsation frequency. The range of the masking was: f_mean ± f_RMS, which is shown in the bottom row of Fig. 6 (the frequencies within the dashed vertical lines). In contrast, the frequency range that was outside the masked region just described was used to recover the signal that was independent of the primary frequency of the pulsations. This signal is shown in the bottom row of Fig. 5.

In the fourth row of Fig. 5, the pulsations for the 1.5 M_⊙ AGB are stable, despite a sudden decrease in the pulsation amplitude just before the 10 year mark. For the less massive 1 M_⊙ AGB with lower surface gravity, in the middle column, the pulsations are still evident, albeit more irregular. There are fluctuations, phase and amplitude shifts in its pulsations, and this is also reflected in its corresponding FFT power spectrum in Fig. 6 – showing many high amplitude neighbouring peaks around a main signal. For the RSG model in the last column, a strong decrease in the pulsation amplitude can be observed in the middle of the simulation. The amplitude stayed low for about two pulsation cycles, then recovering its previous amplitude.

The masking to remove the pulsations enables the observation of the radial component of turbulence and convection within the star, shown in the bottom row of Fig. 5. By comparing the relative intensity of the convective signals, the depths of the zone with significant convective energy flux are shown to be different between the models. The density contrast in the core and atmosphere are different, and, as explained in Sect. 2.1, the inner boundary conditions for different models might vary. Therefore, it is difficult to disentangle and compare any properties of the convective flows in models of different masses. All the representative models show strong overlapping regions where the pulsation and convective signals are strong, indicating that the interplay between pulsations and convection is definitely possible.

Finally, we note that the FFT power spectrum of the 8 M_⊙ model shows a detectable signal just below twice the main pulsation frequency. We verified the signal to be due to an overtone mode in the pulsations. Showing a clear node around 80% of R_* in the power spectrum, the model has the strongest imprint of an overtone in this work. A similar but ambiguous signal is present in the power spectrum of the 1.5 M_⊙ model. However, we attribute this signal to be due to the first harmonic of the main pulsation frequency due to the non-sinusoidal nature of the radius variations. Shock waves in the atmosphere contribute to this signal above the radius R_*. We speculate that, below the radius R_*, interference by convection produces a noticeable signal around twice the pulsation frequency. The relatively large pulsation amplitude makes this feasible. Overall, after investigating the models on an individual basis, strong signals suggesting any overtones are rare. With exception to the 8 M_⊙ model, we favour the interpretation of the second harmonic being the contributor to any detectable signals about twice the pulsation frequency, which were observed in a few of the models.

3.3 Pulsation properties and correlations

The extracted pulsation periods of the stars are plotted against their fundamental stellar parameters in Fig. 7. The plots represent an overview of the current models that we have analysed and the range of stellar parameters we were able to cover so far. Two important trends should be pointed out: (1) increasing pulsation period with radius and; (2) decreasing pulsation period with higher surface gravity. The trends are primarily linked to how the stellar parameters affect the dynamical timescale of the pulsations, and will be presented in further detail alongside other considerations in the following sections.

Fig. 5 Pulsation properties of three representative models. Columns from left to right are outputs of selected 1.5 M⊙, 1 M⊙, and 8 M⊙ models (st28gm05n023, st28gm06n039, and st34gm02n001). The rows, from top to bottom (all showing 15 yr of simulated time): bolometric luminosity; spherically averaged radial velocity; its average in a thin layer beneath the stellar radius; recovered velocity field after masking for the pulsation signal; recovered velocity field after removing the pulsation signal. The solid and dashed lines respectively mark the radius and 80% of the radius of the star. The limits on the colour scales in the second, fourth, and fifth rows are around ±8, ±4, and ±2 [km s−1], respectively, with the colour bars corresponding to the one shown in the top left panel of Fig. 2.

Fig. 6 FFT spectra of the representative 1.5 M⊙, 1 M⊙, and 8 M⊙ models shown in Fig. 5 (st28gm05n023, st28gm06n039, and st34gm02n001) respectively. The FFT power spectra are presented as a function of radial coordinates and frequency (scaled to show the maximum power in between 80–100% of the radius), and with the signal averaged and normalised between the thin layer beneath the stellar radius.

Fig. 7 Derived pulsation periods plotted against selected stellar parameters of the models. Five panels show the pulsation period against (top to bottom reading left to right): surface gravity; stellar radius; averaged effective temperature; stellar luminosity and stellar mass. The last panel illustrates how the models in this work fit into the H-R diagram.

3.3.1 Ritter’s period-mean density relation

Ritter (1879) first established the relationship between the pulsation period and the mean-density of a pulsating homogeneous sphere experiencing adiabatic radial pulsations. The final relation is $P_{\text{puls}}\propto {\overline{\rho }}_{*}^{-0.5}$ and in application to stars, the contraction and expansion of the star would be accompanied by a periodic increase and decrease in effective temperature and luminosity. However, the idea of pulsating stars did not gain much attention until Shapley (1914) suggested radial pulsations, instead of binarity, could explain the variability in Cepheids. Eddington (1918, 1919) then provided further proof and a theoretical framework for such pulsating stars and came to the same conclusions as Ritter (1879); for radial pulsators where the restoring force is the pressure, the pulsation period of a radial mode can be estimated by the time it takes for a sound wave to travel back and forth from the centre to the surface of the star – defined as the sound-crossing time. In the top panel of Fig. 8, we find a good agreement between our models and Ritter’s period-mean density relation. One can derive the pulsation constant, Q, from the period mean-density relation, provided in Fox & Wood (1982): $$P_{\text{puls}}=Q{\left(\frac{{\overline{\rho }}_{*}}{{\overline{\rho }}_{\odot }}\right)}^{-1/2}\to Q=P_{\text{puls}}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{1/2}{\left(\frac{R_{*}}{R_{\odot }}\right)}^{-3/2}.$$(4)

Applying a linear regression on the $P_{\text{puls}}-{\overline{\rho }}_{*}$ plot in Fig. 8, but with ${\overline{\rho }}_{*}\text{/}{\overline{\rho }}_{\odot }$ on the x-axis (using ${\overline{\rho }}_{\odot }=1.3963$ g cm⁻³), provides the coefficients: $$\begin{array}{lll}\log \left(P_{\text{puls}}\right)\hfill & =\hfill & -0.5625\left(\pm 0.0331\right)\log \left({\overline{\rho }}_{*}/{\overline{\rho }}_{\odot }\right)\hfill \\ \hfill & \hfill & +-1.6275\left(\pm 0.2487\right).\hfill \end{array}$$(5)

Therefore, a slightly steeper gradient was obtained against the value derived by Ritter’s relation, as can be seen on the top panel of Fig. 8. The intercept of the derived relation gives the fundamental radial pulsation period of the Sun as 0.0236 ± 0.0182 days, which within its uncertainty agrees with the literature value of 0.033 days (Handler 2013). If a fixed slope of −0.5 was used instead, the corresponding intercept would be 0.0695 ± 0.0111 days, thus significantly overestimating the value of the fundamental radial pulsation period of the Sun. On this premise, the steeper gradient derived is justifiable for the models presented here. However, it is noted that our models follow the $P_{\text{puls}}\propto {\overline{\rho }}_{*}^{-0.5}$ relation well, which is noteworthy given that the models do not fit the underlying assumptions of Ritter’s relation; derived for a homogeneous and spherically symmetric atmosphere experiencing adiabatic pulsations.

Given the dependence of the pulsation constant on the stellar parameters of $Q\propto R_{*}^{-1}{\left(M_{*}\text{/}{R}_{*}\right)}^{0.5}$ for a given period, it is useful to derive the relationship between Q and the stellar mass-to-radius ratio, f_MR = (M_*/M_⊙)/(R_*/R_⊙), as in the bottom panel of Fig. 8. Under the assumption that this relationship is linear and finding the relevant parameters of its relationship through a linear regression, the mass of a pulsating star can be inferred with (by virtue of substituting the relationship into Eq. (4)): $$M_{*}={\left({10}^{-a}{P}_{\text{puls}}{R}_{*}^{b-1.5}\right)}^{1/\left(b-0.5\right)},$$(6)

where P_puls being in days and the M_* and R_* in solar units. The parameters a and b derived from the linear regression on Q and f_MR are presented in Table 1, where different values of the parameters were calculated depending on the different type of models being considered in the regression: AGB stars, RSGs or all together. The reasoning behind this is apparent from the bottom panel of Fig. 8, which indicates two separated regions for the AGB stars and RSGs. In addition, despite showing some degree of linear dependence, the R-squared (R²) values of the linear regressions on Q and f_MR are relatively poor as seen in Table 1.

However, any practical application of Eqs. (5) and (6) to observed data may be limited. We note that the stellar radius value in the relations represents the entropy minimum radius (see Sect. 2.2), which may differ considerably from radius values derived from observations, depending on the optical depth of the visible outer layers.

Fig. 8 Examining the effect of the stellar mean-density on the pulsation periods. Top panel: extracted fundamental pulsation periods against the mean density. Over-plotted is Ritter’s period-mean density relation (dashed line) fitted to the model points and the dash-dotted line is the unconstrained regressed line based on the model points. Bottom panel: pulsation constant, Q, against the stellar mass-to-radius ratio, fMR, where the dotted line is the result of a linear regression applied to all the models in this plot.

3.3.2 P-L relation

The P-L relation of AGB stars is an area that has been extensively studied and has been constrained by observations. Two independent P-L relations are compared with the 3D models for AGB stars, shown in Fig. 9. The first P-L relation was derived from observed carbon-rich Mira variables in the Large Magellanic Cloud (LMC), where a bolometric correction was applied in the K-band to derive the absolute bolometric luminosity (Whitelock et al. 2009). Andriantsaralaza et al. (2022) used improved constraints on the distances to oxygen-rich Mira variables in the Milky Way by making use of observations from Very Long Baseline Interferometry and masers of the AGB stars. They applied the radiative code DUSTY (Ivezic et al. 1999) to derive the bolometric luminosity of the sources. In Fig. 9, the 3D models are in relatively good agreement with both P-L relations, and within the parameter spread and scatter of the observed carbon and oxygen-rich Miras used to derive the respective P-L relations. There is evidence of the universality (that is one can use results from LMC sources in a different galactic environment) of a P-L relation for Miras (see Alvarez et al. 1997; Whitelock et al. 2008; Andriantsaralaza et al. 2022). Nevertheless, one should be cautious in considering such universality, with counter arguments against it (e.g. Barthès et al. 1999; Urago et al. 2020). The existence of some universality is suggested in Fig. 9. However, compared to the observations and their uncertainties, it is sensible to just conclude that our 3D models of AGB stars are within a realistic scatter of real AGB stars.

The RSGs follow a different P-L relation, at the very least having different zero points in the K-band, though the slope may be in agreement with Miras (e.g. Feast et al. 1980; Wood et al. 1983; Jurcevic et al. 2000; Pierce et al. 2000; Kiss et al. 2006; Yang & Jiang 2012; Soraisam et al. 2018). Existence of super AGB stars with M_bol = −8 (Poelarends et al. 2008) also makes the limit between AGB stars and RSGs non-trivial. If the RSG models were to be included in Fig. 9, they would mostly populate the same x-axis range, but with −6 ≥ M_bol ≥ −8. Previous studies suggest the existence of a universal P-L relation for RSGs; however, the uncertainties involved still allow for a wide range of parameter values for the slope and intercept in this P-L relation (see Fig. 7 in Soraisam et al. 2018). It is also risky to convert observations in the K-band into bolometric values for RSGs, as the colour correction formulations are more complicated and uncertain as the spectral evolution of the RSGs needs to be accounted for (see Beasor & Davies 2018; Soraisam et al. 2018) and the distance estimates need to be good to obtain accurate corrections for extinction. Therefore, the comparison between the bolometric P-L relation for RSGs for our models against observations was not attempted here.

Fig. 9 Extracted fundamental pulsation periods against the bolometric luminosity. The solid line is the P-L relation derived from the observed Carbon-rich Miras (plus-signs) in the LMC Whitelock et al. (2009); whereas the dash-dotted line is the P-L relation derived from Oxygen-rich Miras (x-marks) in the Milky Way (Andriantsaralaza et al. 2022). The filled in grey region is the uncertainties in the empirical P-L relations (which overlap one another). Other symbols in the legend indicate the model stellar mass in M⊙ of our AGB models.

3.4 Interaction between convection and pulsation

3.4.1 Interpretation of the spread in the pulsation period

Fluctuations in observed periods and luminosity at a given phase of a pulsation in cool red giants have been detected; several studies suggest that convection plays an important role in the fluctuations, or that the mechanism behind the fluctuations is stochastic in nature (Eddington & Plakidis 1929; Percy & Colivas 1999; Kiss et al. 2006; Percy & Abachi 2013; Cunha et al. 2020). Furthermore, fluctuations can be expected when pulsations interact with the giant convective cells (Schwarzschild 1975; Antia et al. 1984). If this interaction manifests itself in the spread of the pulsation period, σ_puls, the values of σ_puls derived from the models (see Sect. 2.3) help to parameterise the fluctuations. Presented in Fig. 10, there is a good correlation between P_puls and σ_puls. This correlation is to be expected since a longer P_puls is associated with stars of lower surface gravity. Freytag et al. (2017) conclude the size of surface convection cells to be increasing with decreasing surface gravity in the AGB models. Therefore, larger surface convection cells and a more extended atmosphere, both consequently due to a lower surface gravity, could directly cause a larger spread in the pulsation period.

As a result of the positive correlation between P_puls and σ_puls, the relations between σ_puls and the fundamental stellar parameters follow closely with how P_puls depends on the stellar parameters as in Fig. 7. However, one important consideration is how σ∙puis would depend on the number of pulsations observed in the simulation, N_puls. For the period determination process to be reliable, a sufficiently high number of pulsations need to be considered so that the σ_puls value converges to a final value. As seen in Fig. 5, the 1 M_⊙ model had multiple strong peaks in its FFT spectrum, resulting in a larger σ_puls. Values of σ_puls do not necessarily decrease with more observed pulsation cycles in the simulations. Primarily this indicates the fluctuations are consistent throughout the simulations. The relative difference in the σ_puls, whether using a relative or adaptive width to do the weighted averaging in the period determination process, was calculated with $\text{Δ}{\sigma }_{\text{puls}}^{\text{norm}}=\left({\sigma }_{\text{puls}}^{\text{Ada}}-{\sigma }_{\text{puls}}^{\text{Rel}}\right)/{\sigma }_{\text{puls}}^{\text{Rel}}$ and was also plotted against P_puls and N_puls. Where ∆σ_puls is large in the bottom row of Fig. 10, it provides further indication of the interplay between convection and pulsations, as the adaptive width method captures signals surrounding the pulsation signal which are missed by the relative width method. For the 1 M_⊙ AGB stars, convection can lead to substantial signals being present in the FFT power spectrum, in part due to the AGB stars having low surface gravity. On the other hand, RSGs with high ∆σ_puls are indicative of lower pulsation amplitude. Signals due the interference of convection on the pulsations become relatively more appreciable hence would lead to higher σ_puls using the adaptive width method. In this way, Fig. 10 captures an important point that in our models, it is clear that convection and pulsations do interact continuously.

Fig. 10 Dependence of the spread in the pulsation period on observed pulsation properties. The panels show the derived spread in the pulsation period (top row) and relative difference in the spread between using the relative or adaptive method (bottom row), against the pulsation period and the number of pulsations observed in the simulations (left and right columns, respectively).

Fig. 11 Standard deviation in the four periods calculated when the sequence was split into quarters against the spread of the period in the whole sequence. The ratio of the averaged period across the split quarters and the full sequence is coloured in. A similar plot was obtained when the full sequence was split into halves.

3.4.2 Regularity or irregularity of the pulsations

For models with low values of σ_puls in Fig. 10 it does not mean that there is no variability in the pulsations being observed. Examining instances of the regions where the pulsation amplitude significantly decreases in Fig. 5, what follows such events can be a slight phase shift in the pulsations and variation in the amplitude. Thus this may contribute to non-zero values of σ_puls. Another way to identify the variability explicitly was to slice the simulation of models into parts, and repeat the pulsation period extraction procedure on the split parts. This was done on models with at least 3000 time steps. The resulting derived periods are represented in Fig. 11, where a good agreement exists between the scatter of the period in the split sequences against the spread of the period in the full sequence.

It is notable that apart from a few extreme cases, the averaged pulsation period across the split sequences agree with the period derived from the full sequence to about 5%. In addition, it was verified that almost all of the pulsation periods from the split sequences are within the main spread in the period of the full sequence, and when outside, their respective period spreads overlap with each other. The lack of outliers primarily indicates that the pulsation extraction process is reliable, considering the degree of irregularity in the dynamics of the pulsating star. It also shows good regularity of the pulsations despite the variability. Where there are large values in the spread of the period in any of the split sequences, it suggests an event that affects the pulsations. An example of this is the sudden decreases in the pulsation amplitudes observed in the representative models seen in Fig. 5. However, this again highlights the interaction between convection and pulsations, where in certain intervals of the simulation, fluctuations and changes in the period were detected.

Fig. 12 Box plots investigating the choice of stellar radius. Top panel: box plot of the ratio of the pulsation period and sound-crossing time. Bottom panel: box plot of the ratio of the stellar radius and the radius below which the layers are largely in phase. The sound-crossing times and stellar radii are from the six radii definitions tested in the models presented in this work as outlined in Sect. 2.2. Each box plot represents the distribution by drawing the box from first to third quartiles (gi and Q3, respectively) with a horizontal line drawn in between to denote the median, and the height of the box is the interquartile range (IQR). The bottom and top whiskers are bounded by bars denoting the minimum and maximum values, defined as the lowest and highest values within the range min ≥ 1.5IQR − Q1 and max ≤ 1.5IQR + Q3, respectively. Outliers are identified as values being outside of the whisker bounds, denoted by the rhombus symbols in the plot.

3.5 Stellar radius parameter

3.5.1 Entropy minimum as the radius definition

In the analyses presented in this work, choosing a consistent definition for the radius was extremely important, as several physical parameters and properties were derived from it. The multiple chosen definitions of the radius were introduced in Sect. 2.2. Barring the radius estimated from the phase of the pulsation, a sound-crossing timescale was computed for each radius definition. This sound-crossing time was defined as the time it takes for a sound wave to travel back and forth from the centre to the defined radius of the star. We assume the sound-crossing time to be a good first-order estimation of the pulsation period, and first verified that the two timescales correlate linearly to each other. We can equivalently expect the free-fall timescale (over a substantial part of the radius of the star) to be in reasonable agreement with both the sound-crossing and pulsation timescales; by considering the virial theorem (Christensen-Dalsgaard 1993).

The discrepancies in the pulsation periods derived from using the different radius definitions (since the location of the thin layer being considered is different in the weighted averaging, see Sect. 2) were verified to be minimal. The ratio between the pulsation period and the computed sound-crossing times from the different definitions are presented in the top panel of Fig. 12 in the form of box plots. While we expect the sound-crossing times to be a good approximation for the pulsation periods, we note the median values of the box plots in the top panel of Fig. 12 are surprisingly close to 1 − and not with a systematic offset from 1 as initially expected. Ideally, the box plots should have a thin interquartile range with minimal outliers. The entropy minimum definition best matches this case. Using the average temperature of 8000 K and Rosseland optical depths of 10 underestimate the extent of the radius, whereas using Rosseland optical depths of 2/3 and 1 overestimate the radius and the respective box plots have relatively large interquartile ranges. The radius definition used in Freytag et al. (2017) fared well in this test; however, there were outliers which can be attributed to the spatial fluctuations in density and opacity in the stellar atmosphere, which can be mistaken as the boundary for the stellar radius.

To further investigate these statements, box plots of the ratio between the radius derived from the different radius definitions and the averaged radius where the pulsations are in phase are presented in the bottom panel of Fig. 12. We find qualitative agreement and draw the same conclusions with the box plots in the top panel of Fig. 12. Additionally, linear regressions were applied on the relationships being analysed in Fig. 12; the ideal case for the linear regression would simply be (for Y = a + bX) a = 0 and b = 1, with R² = 1. The entropy minimum definition best matched this case in both relationships. These results reaffirm that the radius, as defined by the entropy minimum, was the most consistent and reliable method to calculate the radius in the models. Currently there is no equivalent definition available from an observational stand point. In terms of pulsations, however, we were interested in capturing the outer boundary containing most of the dynamical action. The radius as determined by the entropy minimum definition appears to be the best option in our models.

3.5.2 Radius amplitude due to pulsations

The standard deviation of the radius variation in the entropy minimum definition of the radius, displayed in the right panel of Fig. 1, was calculated and taken as the pulsation amplitude, ΔR_puls. The corresponding ratio ΔR_puls/R_*, hereby called the relative amplitude, is plotted against various parameters in Fig. 13. In general, the relative amplitude follows the trends set out by previous analyses and seems to react strongly to the surface gravity. Against P_puls, the relative amplitude increases with longer pulsation periods (associated with low surface gravity). Following this, the relative amplitude increases with decreasing values of surface gravity, mean-density, and mass-to-radius ratio, specifically for each mass class of the models. It is also worth pointing out that despite low relative amplitudes, the RSGs still produce shocks and strong outward flows in the outer atmosphere.

The relative amplitude increases almost linearly with the quantity v_surfP_puls/R_*, where v_surf is the surface velocity, calculated as the standard deviation value (over time) of the spherically averaged velocity at the radius. This quantifies how well the model represents a perfectly spherically symmetric model with strong radial pulsations. In the bottom right panel of Fig. 13, the three AGB stars that deviate from the linear relation are of low surface gravity with very extended atmospheres. The pulsations, together with the low surface gravity and extended atmosphere, enhance the likelihood for the atmosphere to deviate further away from spherical symmetry.

4 Discussion

4.1 Pulsation constant parameter

Values of Q may help to characterise different pulsation modes (for example, radial or non-radial, or, fundamental or overtone modes) existing in a star. Given its dependence on the stellar parameters shown in Eq. (4) and for a given period, different values of Q indirectly differentiate the stellar interior, especially with regards to the temperature stratification in the star. As suggested in Kiss et al. (2006) however it could be useful to use Q to verify or compare stellar parameters as a check between theoretical models and observations.

Works to derive theoretical values of Q have been produced through usage of 1D models. Xiong & Deng (2007) derived the values of Q analytically from pulsations prescribed by non-local time-dependent convection in 1D for AGB stars (M_* ≤ 3 M_⊙). From their models, the Q values ranged between 0.06 and 0.1 days, with an average of 0.08 and a 1σ scatter of 0.077 ± 0.01 days. In comparison, the 1σ scatters for our models of AGB stars, RSGs and as a whole were respectively calculated to be: Q^AGB = 0.0724 ± 0.0049, Q^RSG = 0.0644 ± 0.0105 and Q^all = 0.0699 ± 0.008 days. Xiong & Deng (2007) had a considerable amount of scatter in the values of Q (see Fig. 1 of Takeuti et al. 2013), but there was a tight correlation between the Q values in the fundamental and overtone modes of their models. Takeuti et al. (2013) makes use of the tight correlation to constrain masses of AGB stars with observed radii, thus presenting a new method to derive masses of observed AGB stars. For this method to be reliable, however, it requires a relatively certain measure of the observed radii and that the theoretical AGB models pulsate in more than one mode (in the fundamental and at least one overtone mode). The lack of overtones in our AGB models means it is not currently possible to apply a similar method and compare with the estimates for the masses of the AGB stars in Takeuti etal. (2013).

Cruzalèbes et al. (2013) derive values of Q for a number of late-type giants. They used observations from the VLTI/AMBER facility (in the K-band) and the MARCS stellar atmosphere models (Gustafsson et al. 2008) to derive new values of angular diameters with higher certainty than in previous studies. The initial input stellar parameters for the MARCS models were effective temperature, surface gravity and stellar mass adopted from the spectral type of the stars in their sample. Due to the multiple options of using different spectral types to infer the input stellar parameters, they confirmed that the derived angular diameters were not significantly affected by the choice of input stellar parameters. Finally, the stellar mass of the targets in their sample was estimated by correlating the temperature and luminosity acquired from the MARCS models to evolutionary tracks in the H-R diagram. From their sample, there were only 4 Mira variables in the fundamental pulsation mode which could be compared to the AGB stars in the present work, which are over-plotted in Fig. 14.

For RSGs, Fadeyev (2012) uses the solution of radiation-hydrodynamics and turbulent time-dependent convection equations to study non-linear oscillations in RSGs with the zero age main sequence (ZAMS) masses in the range of 8 M_⊙ ≤ M_ZAMS ≤ 20_⊙. From their set of models, the relation of their corresponding plot of Q to the stellar mass-to-radius ratio provides the relation log Q = −2.288 − 0.778 log f_MR. We derive the same relation in from our RSG models in this work (M_* ≥ 5 M_⊙) being log Q = −(2.238 ± 0.177) − (0.527 ± 0.094) log f_MR, where the log Q − log f_MR relation is reflected in the bottom row of Fig. 8. Fadeyev (2012) further extended the use of their derived Q to estimate the stellar masses of 7 galactic RSGs with available periods. They made use of the radii determined from MARCS models in Levesque et al. (2005) and Josselin & Plez (2007) and acquired the corresponding stellar masses using their values of Q. From their derived stellar parameters, the corresponding ß values (for RSGs with pulsation periods higher than 100 days) were calculated in order to be compared to our models, and are over-plotted accordingly in Fig. 14.

While there is some overlap, Fig. 14 overall shows there is discrepancy between the estimates of Q of observed stars provided by Cruzalèbes et al. (2013) and Fadeyev (2012), and the models in this present work. Different theoretical models being compared, alongside observational uncertainties involved, contribute to the discrepancies. Furthermore, factoring in the high uncertainties in the radius estimations inferred from observations, as well as further differences in the definition of the radius of a star (see Sects. 3.5.1 and 4.4), the discrepancies between our work here and the values derived from observations in Cruzalèbes et al. (2013) and Fadeyev (2012) are to an extent unavoidable. Both of the works used the Rosseland radius definition (τ_Ross = 1) in the MARCS models to define the radii being used. However, any radii calculated via the MARCS models should be taken with a level of caution, as MARCS assumes static atmospheres, whereas the stars being compared here are dynamical in nature, evident from the pulsations. Hence it is challenging to compare the works here, especially as the parameters were indeed derived from different methods, which involves their own uncertainties; however, Fig. 14 represents an exploratory attempt to find comparisons in literature and observations, and show that the pulsation constants of our models are within the general range derived from available observations.

Fig. 13 Effect of the pulsations on the stellar radius. The panels show the relative amplitude due to pulsations against (from top to bottom reading left to right): surface gravity; effective temperature; mean density; mass-to-radius ratio; pulsation period and the surface velocity multiplied by the pulsation period over the stellar radius.

Fig. 14 Logarithm of the pulsation constant, Q, against the logarithm of pulsation period. To compare against the models, 4 Mira variables (rhombuses) and 7 RSGs (pentagons) are included from Cruzalèbes et al. (2013) and Fadeyev (2012), respectively. Uncertainty bars have been excluded for the models to preserve the clarity of the plots.

4.2 Difference in models

Section 2.1 introduced the general setup of the models, further developed from previous works in 3D simulations of pulsating AGB stars (see Freytag & Höfner 2008; Freytag et al. 2017). As part of the work presented here, pulsations extracted from different generations of the CO5BOLD code were investigated. Higher resolution models introduced additional finer convective structures (Freytag et al. 2017) and slightly more efficient convection in the models, but the period of the pulsations was consistent with lower resolution models. Temporally, longer running simulations help to investigate the fluctuations in the pulsations more reliably – and smaller sampling time steps do not reveal anything new.

The inclusion of radiative pressure in the stellar interior becomes more important with higher stellar mass and inner temperatures (see Sect. 2.1 and Goldberg et al. 2022). After examining the radiation-to-gas pressure ratio in the interior, the missing radiative pressure is expected to affect the deep stellar interior of RSGs, and could enhance convectively unstable regions in the core of the RSG models. In Figs. 8 and 14, there is a 6M_⊙ RSG model that seemingly identifies as an outlier when compared to other RSGs. This model (st36g00n04) was one of the first attempts to produce a representative model of Alpha Orionis (Betelgeuse). The model has realistic temperature and surface gravity but smaller mass and radius than the inferred stellar parameters of Betelgeuse. Temperature and gravity are at the upper end of our model set (see Fig. 7.), while the numerical resolution is comparatively low. This leads to relatively inefficient convection not carrying energy up from the core region, where radiative energy transport dominates, but only in the outer layers of the star. The entropy profile formally indicates convective stability of the inner region, while overshooting convective velocities still can reach very deep down. In this case, we believe that the inclusion of radiation pressure would deepen the convection zone. The relatively shallow convection zone might have consequences for the excitation of the pulsations. But it should have only a minor impact on the pulsation period, because for the period the outer layers, where the sound speed is low, are of dominating influence. Thus, st36g00n04 serves as an extreme case to highlight possible limitations in our current set of RSG models. However, it was verified that this model, too, is pulsating in the fundamental mode. Hence, its results still appear valid and useful to maximise the range in our parameter space in the analyses presented in this work.

The current generation of models have only explored a relatively tight region of parameter space in the mass-to-radius ratio, f_MR (see Sect. 3.3.1), in each mass class, indicated in the bottom panel of Fig. 8. The reasoning for this is different for the AGB and RSG models. For AGB models, the choice of stellar parameters was partially guided by DARWIN models that produce stellar winds driven by radiation pressure on dust (see Eriksson et al. 2014; Bladh et al. 2019a,b). As the escape velocity from a star scales with v_esc ∝ (M_*/R_*)^0.5, a lower f_MR was preferred to generate stellar winds easily. This is since matter can be lifted up to distances from the star, where low temperatures permit dust condensation. In global RSG models, consideration towards the size of the convection cells at the surface and the ability to resolve them becomes important; where the size of granules, or convection cells, is proportional to the pressure scale height (see for instance Freytag et al. 1997, and references therein). Larger pressure scale heights are then preferred since this would allow surface features to be resolvable. In relation to the ratio of the radius and pressure scale height at the photosphere: R_*/H, with $H\propto 1/g\propto R_{*}^2/M_{*}$, the ratio becomes R_*/H ∝ M_*/R_*. Therefore, the preference of larger pressure scale heights has limited the ranges of the f_MR in simulations of the RSGs. In general however, going towards higher ratios of f_MR requires higher spatial resolution, which is more computationally expensive.

4.3 Convection and pulsation interaction

Section 3.4 highlighted the evidence of the interaction between convection and pulsations in the models. This interaction has several effects on the atmosphere, importantly creating inhomogeneities on the surface of the photospheric layer.

The surface inhomogeneities due to convection may interact with the pulsations and atmospheric shock waves. This interaction produces a network of shock waves in the upper atmosphere (Freytag et al. 2017). In this sense, the pulsations and shock waves can carry forwards any spatial inhomogeneities on the surface, creating patchy and irregular outflows. Thus it is thought the size of surface convective cells and how the pulsations interact with convection are significant factors for the production of irregular and patchy outflows. The interaction also causes variability in the pulsation period and amplitude.

While the study of the interaction was not presented here, it is evident decoupling the pulsations and convection is difficult. The reason for the difficulty in investigating the non-linear interplay between convection and pulsations is that there is no explicit parameterisation for how the convection affects the pulsation, and vice versa, in the CO5BOLD code. Convective signals still remain substantial despite applying spherical averages to obtain the radial velocity: observed in the first and third rows of Fig. 5. In fact, the pulsations heavily modulate the convection energy flux (the sum of the enthalpy and kinetic energy fluxes) at every pulsation cycle. Reacting against this, the stratification in the stellar interior and atmosphere may experience perturbations in the form of pressure fluctuations. The mass, momentum and energy fluxes are affected, and the degree of this disturbance can influence the global behaviour of the pulsations. We interpret sudden amplitude change in the pulsations, as observed for the 1.5 M_⊙ and 8 M_⊙ models in the top row Fig. 2, to be due to interference by convection. Given the high Mach numbers (velocities exceeding 1 km s⁻¹) of the convective flows and modulation due to pulsations, quantifying the interaction between convection and pulsations remains challenging.

Figures 10 and 11 investigate the relationship between fluctuations in the pulsation period and stellar parameters. The investigations represent a phenomenological way to study the effect of convection on the pulsation period. For models that experience regular, high amplitude pulsations (interpreted as models with low σ_puls), convection still causes some variability in the period following events where pulsation amplitudes suddenly decrease. Capturing the variations and asymmetries due to the non-linear interaction between convection and pulsations is important in reproducing more realistic dynamical models and the mass loss mechanism (see Liljegren et al. 2016).

4.4 Stellar radius definition

Works determining the radius of AGB stars and RSGs from observations fundamentally rely on how the radius, that is the boundary of the star, is first defined. Historically there have been competing ideas on which definitions work best, especially for stars with cool extended envelopes (Haniff et al. 1995). Furthermore, the extent of the surface of cool red giants carries high uncertainty due to time-variability, non-sphericity, and extensions of the surface of the star arising from the dynamical atmosphere (Arroyo-Torres et al. 2015). Therefore, it becomes non-trivial to decide how to define the radius from observations and how to compare with the radii obtained in theoretical models.

It is difficult to compare the radii of AGB stars and RSGs obtained from the CO5BOLD simulations against radii estimated from observations. Fundamentally, this comes down to the radius definition being used as suggested in Sect. 3.5.1. From the observational view, the dynamical nature of the atmosphere and the molecular and dust formation shells surrounding cool evolved stars (e.g. Wittkowski et al. 2007; Khouri et al. 2016; Ohnaka et al. 2017; Höfner & Olofsson 2018) evidently present themselves as uncertainties to the true radius of the star. Analogous to the different conditions which currently exist to describe the radius of the star in the CO5BOLD models, radii from the observational stand point also have various definitions, being mostly dependent on the wavelength observed and what is possible technologically. The distance measurement to target sources is also a major uncertainty in the observed radii of cool red giants (see Baschek et al. 1991; Haniff et al. 1995; Chiavassa et al. 2020). Despite the complications, it is however reassuring that for the AGB models, Fig. 9 shows that the over-plotted P-L relations are in good agreement with the scatter of the observed Miras, indicating that the stellar parameters in, or acquired from, the models are in the region of those of real stars.

Section 3.5 investigated which definition for the radius yielded the most reliable radius of the modelled star. It is important to state that our primary concern was to decide on a definition that would capture the shell where most of the dynamic action (contained within the star), particularly pulsations, took place. The entropy minimum definition, marking the region where the atmosphere becomes convectively stable, appeared to be the best definition in Fig 12. Several different stellar parameters, such as the effective temperature, mean-density and mass-to-radius ratio, can be derived from the radius definitions, hence finding the most reliable definition is crucial. Using the entropy minimum definition also yielded the closest agreement to the $P_{\text{puls}}-{\overline{\rho }}_{*}^{-0.5}$ relation as in Fig. 8. This reassures that the radius definition we adopted is the best option amongst the other definitions presented in this work. The closest equivalent observational definition of the radius to this would be the very inner boundary of the photometric radius of the observed star. This requires either a combination of high precision interferometry and distance estimates to the star itself, or precise radiative transfer modelling of the luminosity and temperature of the star to obtain estimates for the photometric radius of the star.

5 Summary and conclusions

Global 3D radiation-hydrodynamics models of AGB and RSG stars are an essential tool for studying their pulsation properties and the excitation mechanism. The 3D models computed with the CO5BOLD code were used to calculate the fundamental pulsation frequency and other critical pulsation properties. A total of 75 models was analysed in this work, comprising several generations of models. The convection and pulsations are both emergent properties of the radiation-hydrodynamics simulations. This means that no parametric treatment of convective energy transport is necessary, in contrast to classical 1D pulsation models. Despite the non-linear nature of the interplay between convection and pulsations, in all of the models radial pulsations can clearly be observed and the pulsation periods in the fundamental mode were extracted. A good correlation was found for the $P_{\text{puls}}-{\overline{\rho }}_{*}^{-0.5}$ relation in the models, confirming that the pulsations are in the fundamental radial mode. This agrees with current theory that for LPVs, the variability is attributed to low-order large-amplitude pulsations, existing as standing waves.

In their dependence on fundamental stellar parameters, the pulsation periods agree with expectation, such as a lower surface gravity resulting in longer periods. The AGB models agree well with available observations in the P-L relation, indicating that the stellar parameters and the pulsation properties are realistic.

Applying regression models led to the mass estimate formula in Eq. (6), as a function of the fundamental pulsation period and stellar radius. The radius value to be used in Eq. (6) must be consistent with the radius used when the formula was derived. We found that identifying the innermost local entropy minimum in the atmosphere provided the best measure for the stellar radius in the context of pulsation analysis. The purpose of identifying the best radius definition was to explicitly capture the boundary of dynamics in the stellar interior, thus where the pulsations play a dominant role.

It is clear from the simulations that convection actively affects the pulsation and that their interplay causes period variability. There is evidence that a reverse effect exists, in the sense that the pulsations also impact the convection. The deep convective regions and large convection cells characteristic of evolved red giants make it important to understand the interplay with pulsation, and ultimately how the interplay affects the massive stellar winds and the mass loss process. Further and more in-depth studies on the interaction of convection and pulsations in the 3D models remain on the agenda for future work. It proves to be a difficult topic to investigate, since the effect of one cannot be isolated, or decoupled, from the other.

Acknowledgements

This work is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 883867, project EXWINGS) and the Swedish Research Council (Vetenskapsrådet, grant number 2019-04059). Most of the computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX and NSC in Sweden, supported by the Swedish Research Council (Vetenskapsrådet). Some simulations were performed on resources provided by PSMN in Lyon, France, and MPA in Garching, Germany.

Appendix A Model parameters and results

References

All Tables

All Figures

Fig. 1 Entropy profile of the standard model st26gm07n001. Left panel: temporal mean value of the entropy as a function of the radial coordinates, with the shaded area being the corresponding standard deviation and the zoomed-in axes shows the first local entropy minimum after the first instance of a spherically averaged temperature of 8000 K away from the centre of the star. Right panel: count values of where the entropy minimum was located over a simulation (in bins of size 2 R⊙).

In the text

Fig. 2 Diagnostics of the spherically averaged radial velocity for the standard model st26gm07n001. Top left: radial velocity as a function of time and radial coordinates, where the solid and dashed lines mark the radius and 80% of the radius, respectively. Bottom left: average and standard deviation of the radial velocity (in blue and shaded in grey, respectively) within the layer of 80–100% of the radius. Top right: FFT power spectrum (values are clipped to the maximum power calculated in between the 80–100% of the radius) as a function of the frequency and radial coordinates. Bottom right: averaged and normalised FFT power spectrum across the aforementioned layer. Here the shaded region is the extent over which the weighted averaging (see text) to get a characteristic pulsation frequency was done, while the solid and dashed black vertical lines indicate the computed characteristic frequency and its spread, respectively.

In the text

Fig. 3 Radial velocity profiles for the standard model st26gm07n001. Left: snapshots of the spherically averaged radial velocity as a function of the radial distance, covering 4 yr of simulated time (about three pulsation cycles). Right: radial velocity computed at each cell of the computational box determined from a snapshot at time t = 14.2 [yr], presented to show the number of points (N) having the radial velocity at a certain distance from the centre of the box. The vertical lines mark the location of the stellar radius.

In the text

	Fig. 4 Comparison between the different width methods, relative (Rel) or adaptive (Ada) in superscript, in the derived pulsation periods and spreads (left and right panels, respectively). Two models that have extremely extended atmospheres were omitted from the right panel to preserve clarity of the plot (as the values are around 0.8 [yr]).
In the text

Fig. 5 Pulsation properties of three representative models. Columns from left to right are outputs of selected 1.5 M⊙, 1 M⊙, and 8 M⊙ models (st28gm05n023, st28gm06n039, and st34gm02n001). The rows, from top to bottom (all showing 15 yr of simulated time): bolometric luminosity; spherically averaged radial velocity; its average in a thin layer beneath the stellar radius; recovered velocity field after masking for the pulsation signal; recovered velocity field after removing the pulsation signal. The solid and dashed lines respectively mark the radius and 80% of the radius of the star. The limits on the colour scales in the second, fourth, and fifth rows are around ±8, ±4, and ±2 [km s−1], respectively, with the colour bars corresponding to the one shown in the top left panel of Fig. 2.

In the text

	Fig. 6 FFT spectra of the representative 1.5 M⊙, 1 M⊙, and 8 M⊙ models shown in Fig. 5 (st28gm05n023, st28gm06n039, and st34gm02n001) respectively. The FFT power spectra are presented as a function of radial coordinates and frequency (scaled to show the maximum power in between 80–100% of the radius), and with the signal averaged and normalised between the thin layer beneath the stellar radius.
In the text

	Fig. 7 Derived pulsation periods plotted against selected stellar parameters of the models. Five panels show the pulsation period against (top to bottom reading left to right): surface gravity; stellar radius; averaged effective temperature; stellar luminosity and stellar mass. The last panel illustrates how the models in this work fit into the H-R diagram.
In the text

Fig. 8 Examining the effect of the stellar mean-density on the pulsation periods. Top panel: extracted fundamental pulsation periods against the mean density. Over-plotted is Ritter’s period-mean density relation (dashed line) fitted to the model points and the dash-dotted line is the unconstrained regressed line based on the model points. Bottom panel: pulsation constant, Q, against the stellar mass-to-radius ratio, fMR, where the dotted line is the result of a linear regression applied to all the models in this plot.

In the text

Fig. 9 Extracted fundamental pulsation periods against the bolometric luminosity. The solid line is the P-L relation derived from the observed Carbon-rich Miras (plus-signs) in the LMC Whitelock et al. (2009); whereas the dash-dotted line is the P-L relation derived from Oxygen-rich Miras (x-marks) in the Milky Way (Andriantsaralaza et al. 2022). The filled in grey region is the uncertainties in the empirical P-L relations (which overlap one another). Other symbols in the legend indicate the model stellar mass in M⊙ of our AGB models.

In the text

	Fig. 10 Dependence of the spread in the pulsation period on observed pulsation properties. The panels show the derived spread in the pulsation period (top row) and relative difference in the spread between using the relative or adaptive method (bottom row), against the pulsation period and the number of pulsations observed in the simulations (left and right columns, respectively).
In the text

	Fig. 11 Standard deviation in the four periods calculated when the sequence was split into quarters against the spread of the period in the whole sequence. The ratio of the averaged period across the split quarters and the full sequence is coloured in. A similar plot was obtained when the full sequence was split into halves.
In the text

Fig. 12 Box plots investigating the choice of stellar radius. Top panel: box plot of the ratio of the pulsation period and sound-crossing time. Bottom panel: box plot of the ratio of the stellar radius and the radius below which the layers are largely in phase. The sound-crossing times and stellar radii are from the six radii definitions tested in the models presented in this work as outlined in Sect. 2.2. Each box plot represents the distribution by drawing the box from first to third quartiles (gi and Q3, respectively) with a horizontal line drawn in between to denote the median, and the height of the box is the interquartile range (IQR). The bottom and top whiskers are bounded by bars denoting the minimum and maximum values, defined as the lowest and highest values within the range min ≥ 1.5IQR − Q1 and max ≤ 1.5IQR + Q3, respectively. Outliers are identified as values being outside of the whisker bounds, denoted by the rhombus symbols in the plot.

In the text

	Fig. 13 Effect of the pulsations on the stellar radius. The panels show the relative amplitude due to pulsations against (from top to bottom reading left to right): surface gravity; effective temperature; mean density; mass-to-radius ratio; pulsation period and the surface velocity multiplied by the pulsation period over the stellar radius.
In the text

	Fig. 14 Logarithm of the pulsation constant, Q, against the logarithm of pulsation period. To compare against the models, 4 Mira variables (rhombuses) and 7 RSGs (pentagons) are included from Cruzalèbes et al. (2013) and Fadeyev (2012), respectively. Uncertainty bars have been excluded for the models to preserve the clarity of the plots.
In the text