© The Authors 2025

1. Introduction

Red supergiants (RSGs) mark the end stages of stellar evolution in a certain mass range of massive stars. Some defining features of RSGs are that they are located at the cool edge of the Hertzsprung-Russell diagram (HRD) and have a deep convective and hydrogen-rich envelope that extends from several hundred to thousands of solar radii in size. Observationally, RSGs are known to show (periodic) brightness variations (Kiss et al. 2006; Percy & Khatu 2014) – some well-known variable RSGs are Betelgeuse (α Ori), VY CMa, TV Gem, VX Sgr, and S Per. These brightness variations have been attributed to radial pulsations in the envelope of the RSGs (Stothers 1969; Stothers & Leung 1971; Heger et al. 1997; Wood et al. 1983; Guo & Li 2002; Yoon & Cantiello 2010) and are thought to be driven by the κ mechanism. Alternative explanations for the irregular variability include the brightness variations caused by the large convective cells present in the envelopes of RSGs (Schwarzschild 1975; Antia et al. 1984).

Most stars that end their lives as RSGs explode in a supernova (SN). Supernovae that show a clear signature of hydrogen are classified as Type II SNe (Minkowski 1941) and are connected to explosions of massive stars with a hydrogen-rich envelope (Smartt 2009, 2015). Type II-P SNe are a subclass of Type II SNe whose light curves show a prominent plateau phase of ≈100 d. Type II-L SNe are a second class of hydrogen-rich SNe with a linearly declining light curve (Barbon et al. 1979). Because of the low numbers of Type II-L SN observations (Smartt 2009), it was believed that there is a clear distinction between II-P and II-L SNe (Patat et al. 1993, 1994; Arcavi et al. 2012). However, as the numbers of SN observations increased, there seemed to be a less clear distinction between these two categories and more of a continuum between II-P and II-L SNe, characterized by different decline rates of the light curve during the first ≈100 d (Anderson et al. 2014; Faran et al. 2014a,b; Sanders et al. 2015). Additionally, SNe with a fast-declining light curve are found to be more luminous compared to SNe with lower decline rates Anderson et al. (2024).

Connecting SN properties and the light-curve decline rate to their exact progenitors remains a challenge, though significant progress has been made in recent years. Using archival imaging, it is possible to identify RSGs as the direct progenitor stars of the classical Type II-P SNe that have low decline rates (Smartt 2009). The progenitors of SNe with a high decline rate, i.e., classically referred to as Type II-L SNe, are not well known due to their scarcity. There seems to be a trend that yellow supergiants might be possible progenitors (Van Dyk 2016). However, it should be noted that there are exceptions to this trend, and Valenti et al. (2016) claim that there is no clear distinction between the progenitors of SNe with faster or slower decline rates.

Numerical models computing the SN light curves from pre-SN stellar structures show that the explosions of RSGs are indeed able to reproduce observed light curves of hydrogen-rich SNe with flat light curves (e.g., Litvinova & Nadezhin 1983; Chieffi et al. 2003; Bersten et al. 2011; Dessart et al. 2013; Morozova et al. 2015). Producing fast-declining SNe from numerical models can be achieved by stellar structure with less massive hydrogen-rich envelopes (Grassberg et al. 1971; Young & Branch 1989; Blinnikov & Bartunov 1993; Schlegel 1996), although reducing the envelope mass alone fails to reproduce the observed luminosities (e.g., Hillier & Dessart 2019). Potential sources for reduced envelope masses can be strong mass loss in the RSG phase (Moriya et al. 2011; Georgy 2012; Fuller & Tsuna 2024), mass transfer in binary systems (Podsiadlowski et al. 1992; Nomoto et al. 1995; Dessart et al. 2024), or even eruptive mass loss shortly before the SN (Smith & Arnett 2014; Clayton 2018). The presence of strong mass loss is backed up observationally, as many SNe with fast-declining light curves have a larger peak brightness and often show signs of interaction with circumstellar material (CSM) (Smith 2014; Bostroem et al. 2019; Zhang et al. 2024; Jacobson-Galán et al. 2024).

Supernova light curve calculations require an input stellar structure that usually comes from a one-dimensional (1D) stellar evolution calculation. Generally, hydrostatic equilibrium is assumed in such stellar evolution calculations. This is also the case for pre-SN structures that enter SN light curve calculations¹. Large-scale radial pulsations in RSGs are usually ignored. However, the envelope structure of a pulsating RSG might significantly deviate from a hydrostatic envelope, depending on the amplitudes, but the structure of the hydrogen-rich envelope is critical in determining the shape of the resulting SN light curve. One notable exception is the work by Goldberg et al. (2020), which studies the effect of the pulsations of RSGs just before the onset of core collapse. They analyzed at which frequencies the RSGs would pulsate and then triggered these pulsations just before the SN via velocity perturbations in the envelope. Here, we aim to simulate the pulsations self-consistently whenever the envelopes of RSGs become dynamically unstable.

In this work, we study the radial pulsations in the envelopes of RSGs, their impact on the envelope structure, and the light curves of the resulting SNe. First, we describe the computational tool that we use in Sect. 2. This includes a 1D stellar evolution code to study the pulsations, a dust radiation code to predict light curves of the pulsating RSGs, a semi-analytic code to predict the core-collapse outcome, and a 1D hydrodynamic SN code to compute the SN light curves. In Sect. 3, we analyze in detail the pulsations of RSGs, their origin, their driving mechanism, and their effects on the envelope structure. Then, we model the explosion of the pulsating RSG and compute the resulting SN light curve in Sect. 4. In Sect. 5, we discuss some uncertainties of our models and the implications for SN light curve modeling, before summarizing and concluding in Sect. 6.

2. Methods

In this section, we describe the methods used to compute the pulsating RSG and the SN light curves. The stellar-evolution models are described in Sect. 2.1. Sect. 2.2 summarizes the methods used to model the pulsations of the RSG. The dust computations are explained in Sect.2.3. Lastly, we describe the SN light curve modeling in Sect. 2.4.

2.1. Stellar evolution models

We evolved an initially 15 M_⊙ stellar model from the zero-age main sequence all the way to core collapse, using the Modules for Experiments in Stellar Astrophysics (MESA) stellar-evolution code revision 10398 (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) with the same setup as in Schneider et al. (2021) and Laplace et al. (2025b). In particular, we used a mixing-length parameter α_MLT = 1.8 and employ the Ledoux criterion to model convection. We used MLT++ for better numerical stability because it increases the efficiency of energy transport in convective regions that are close to the Eddington luminosity. Convective-boundary mixing was implemented using step-overshooting of 0.2 pressure-scale heights on top of convective core hydrogen and helium burning. There is no convective-boundary mixing in all other convective shells. Additionally, we used semi-convection with an efficiency factor of α_sc = 0.1.

The model was evolved to core-collapse, defined when iron-core infall velocity exceeds 950 km s⁻¹. At the end of core carbon burning, defined as the moment when the core carbon mass fraction is lower than 10⁻⁴, the star appears as an RSG with a mass of 12.2 M_⊙, a radius of 1024 R_⊙, a luminosity of 1.12 × 10⁵ L_⊙ and an effective temperature of 3304 K. The bottom of the convective envelope is at a mass coordinate of 5.23 M_⊙. This is the starting point for the subsequent dynamical modeling of the RSG.

2.2. Modeling of the pulsations

We used the model of the 15 M_⊙ RSG at the end of core carbon burning as the starting point for the dynamical calculations. First, we enforced hydrostatic equilibrium throughout the entire model. The default settings of MESA for massive stars, as found in inlist_massive_defaults, lift the assumption of hydrostatic equilibrium in layers exceeding a temperature of 10⁸ K when the time step is smaller than 0.1 yr. This way, the contraction of the core during the final burning stages and the onset of core collapse can be directly modeled. The temperature cutoff ensures that the hydrodynamical mode is enabled only in the core and not in the envelope. After forcing the model into hydrostatic equilibrium, we excised the core at a mass coordinate of 4.6 M_⊙. This is about 90% of the helium core mass when defining the helium core by a hydrogen abundance of less than 0.1. The exact mass coordinate at which we remove the core has a negligible effect on the subsequent simulations, as long as it is within the helium core and outside the helium-burning shell (see Appendix B). It is necessary to enforce hydrostatic equilibrium before excising the core to have a static inner boundary condition that does not have any long-term effect on the dynamical calculations.

Once the velocity was forced to zero throughout the star and the core is removed, we enabled the hydrodynamic modeling in MESA throughout the entire simulated domain. We used a basic setup for the hydrodynamic modeling that is similar to the setup used in Clayton et al. (2017) and Bronner et al. (2024). We used MESA’s implicit hydrodynamic capabilities as described in Paxton et al. (2015) that make use of artificial viscosity. For the outer boundary condition, we used the default setting from MESA; that is, modeling an atmosphere using the Eddington T − τ relation (Eddington 1926; Paxton et al. 2013). To ensure that we could resolve the dynamical evolution in the envelope, we limited the time steps to 10⁻² yr. Finally, we turned off any change in abundances from nuclear burning, while still keeping the energy generation rate. This keeps the star in a chemically stable configuration without any nuclear timescale evolution and is justified because of the relatively long timescale of the chemical evolution compared to the dynamical evolution of the envelope.

We ended the simulation after 300 yr, which is much longer than the initial thermal timescale of the envelope of about 11 yr. Longer simulations do not introduce any long-term changes because we disable the conversion of elements via nuclear burning, and we used a static inner boundary condition.

There are several limitations to our pulsation model, most of which originate from the 1D nature of the model and the associated lack of 3D effects on the pulsations. These limitations are presented and discussed in detail in Sect. 5.7. For the models of the pulsating RSG, we in particular keep using MLT++ and a mixing-length parameter of α_MLT = 1.8 for the envelope. Work exists that suggests increasing the value of α_MLT in the envelope of RSGs to produce SNe that align better with observations (e.g., Dessart et al. 2013; Paxton et al. 2018; Goldberg et al. 2022a). This is also addressed in Sect. 5.7.

2.3. Dust modeling

Some RSGs are believed to be surrounded by a dust shell (Verhoelst et al. 2009). Such a dust shell might originate from the mass loss via winds and could be enhanced by pulsations and eruptive mass loss (Bonanos et al. 2024; de Wit et al. 2024). We used the dust radiation-transport code DUSTY (Ivezić & Elitzur 1997, 1999) to compute light curves in different broad-band filters of the RSG model, assuming it is surrounded by a dust shell. For the dust shell, we used common assumptions for RSGs. The dust is composed of warm silicates (Ossenkopf et al. 1992), and the grain-size distribution follows a power law with exponent −3.5 and minimum and maximum grain size of 0.005 μm and 0.25 μm (Mathis et al. 1977). The dust temperature at the inner boundary, r₁, was set to 800 K and the density distribution of the dust shell follows a power law with exponent −2. For the total extent of the dust shell, we use a relative thickness of r_out/r₁ = 1000. The total optical depth, τ, defined at 0.55 μm, is treated as a free parameter and is varied between 0.1 and 100 for each DUSTY model. These assumptions might not be valid for a post-core-helium burning RSG that we consider. Convective mixing after core helium burning could have changed the envelope composition and hence the dust chemistry. Additionally, some RSGs are surrounded by exceptionally large dust grains (Scicluna et al. 2015). Because we are mostly interested in the qualitative effects of a dust shell around a pulsating RSG, we refrain from a more detailed analysis of the dust parameters, but study this in Laplace et al. (2025a), hereafter Paper II.

One main input for the dust radiation calculations is the spectral shape of the illuminating radiation source (i.e., the pulsating RSG). We used the MARCS models (Gustafsson et al. 2008) to describe the spectrum of the RSG. We selected the subset of spherical 5 M_⊙ models² with standard abundances, solar metallicity, and a microturbulence parameter of 5. This subset of MARCS models is available on a well-defined grid of T_eff and log g.

We computed DUSTY+MARCS (DM) models as described above. For each DM model, we computed the absolute V- and K-band magnitudes³ from the dust-attenuated spectral energy distribution. To construct light curves of the RSG, we linearly interpolated between the DM model based on the T_eff and log g of the RSG (for more details on the interpolation, see Appendix D). The DM models only cover temperatures of T_eff > 2400 K, such that we are unable to predict the parts of the light curve where the RSG is cooler than 2400 K. For temperatures cooler than 2400 K, we do not make any predictions, because extrapolation may introduce large uncertainties. The final free parameter of the light curves is the total optical depth, τ₀, at the beginning of the pulsation cycle. This parameter corresponds to the total dust mass around the RSG and cannot be determined a priori. The total optical depth τ₀ can be converted to the total dust mass. This corresponds to dust masses between 10⁻⁴ and 10⁻² M_⊙ for τ₀ = 0.5 − 10, assuming a dust grain density of 3 g cm⁻³. We do not assume any extinction and reddening for these theoretical predictions.

2.4. Supernova light curve models

To determine the SN explosion parameters of the RSG, we employed the semi-analytic SN code by Müller et al. (2016) with the same calibrations as Schneider et al. (2021) and Temaj et al. (2024). Based on our stellar model at the onset of core collapse, applying the code results in an explosion energy of E_exp = 1.60 B, a mass of M_NS = 1.502 M_⊙ for the compact object remnant, and a mass M_Ni = 0.086 M_⊙ of synthesized ⁵⁶Ni. We employed these explosion parameters to compute the SN light curve with the 1D open-source Lagrangian radiation hydrodynamics code SNEC (Morozova et al. 2015) version 1.0.1. SNEC assumes local thermodynamic equilibrium (LTE) for computing the light curves, which is a reasonable assumption during the plateau phase but is less accurate for the shock breakout and nebular phases, which are known to be affected by nonthermal effects. To compute the SN light curves of our progenitor model at different pulsation phases, we mapped each stellar structure on a Lagrangian grid containing 1000 points. We used the thermal bomb model implemented in SNEC to compute the propagation of a shock through the stellar structure. After excising the mass M_NS from the grid, we injected enough energy in the 0.1 M_⊙ above this mass coordinate to reach the desired final energy, E_exp. We simulated strong mixing of nickel in the ejecta by spreading M_Ni of nickel from the inner boundary to 90% of the final mass. We kept the other properties, such as the opacity, equation of state, and ionization states, to the default values for II-P SNe in SNEC and followed the shock propagation until about 200 d after the explosion.

3. Pulsations of red supergiants

In this section, we present the results of the pulsation calculations. First, we show the growth of the pulsations in Sect. 3.1 and describe the pulsation mechanism in Sect. 3.2. Then, we show how the pulsations change the density structure of the RSG in Sect. 3.3, and how such a pulsating RSG might appear observationally by applying a dust shell model around the RSG and computing photometric light curves in Sect. 3.4.

3.1. Growth of pulsations and envelope restructuring

Once we turn on the hydrodynamic mode of MESA at the end of core carbon burning, the RSG naturally starts to pulsate radially (Fig. 1). The amplitude of the radial pulsations grows exponentially for t ≲ 30 yr. We computed the growth rate, η, as defined in Yoon & Cantiello (2010) via η = |v(t₀ + P)/v(t₀)|, where v(t₀) is a relative maximum of the surface velocity and P ≈ 950 d is pulsation period. We find η ≈ 1.87. This is lower than the value of η = 2.15 from Yoon & Cantiello (2010, their equation 1), but still within their observed spread (c.f. Yoon & Cantiello 2010, Fig. 3).

Fig. 1. Radial pulsation of the RSG. The radius evolution is shown in panel a with a zoom-in on one pulsation cycle in panel b. The dash-dotted black line indicates the radius of the hydrostatic model. Panel c shows one pulsation cycle in the HRD. The arrow indicates the direction of the loop in the HRD, and the markers are spaced equally in time every 1/20 of the pulsation period. The hydrostatic model, the time-averaged effective temperature and luminosity, and their uncertainties are shown as well. For comparison, we show the typical uncertainty of the luminosity and effective temperature of Type II SN progenitors from Smartt (2015) as a gray marker. The colored markers labeled A–F in panels b and c indicate characteristic stages during the pulsation at which we calculate SN light curves. The pulsation phase, ϕ, is defined to be zero at maximum expansion.

At a time of around 32.5 yr, there are drastic changes within the entire envelope that cause changes in both the pulsation amplitude and period (i.e., just after peak radius in Figs. 1a and 2a). Figure 2c shows that the density near the bottom of the convective envelope increases. At the same time, the entropy decreases. This change is caused by a catastrophic cooling event at the surface of the extended envelope. Initially, the pulsation amplitude grows exponentially, oscillating around the equilibrium radius. At the same time, the partial ionization zones of hydrogen and helium periodically move inward and outward (Fig. 2b). This means that ionization energy is periodically released (i.e., usually referred to as recombination energy) during the expansion, and energy from the pulsations is used to ionize hydrogen and helium during the contraction. The amount of released recombination energy during each pulsation cycle increases exponentially because the locations of the partial ionization zones, more specifically their mass coordinates, vary with an exponentially growing amplitude. The total released recombination energy per cycle is directly proportional to the difference in the mass coordinate of the partial ionization zones at maximum compression compared to maximum expansion. At t = 32.5 yr, the envelope extends to such large radii that about 2 M_⊙ of the envelope are completely neutral. Additionally, the recombination front of doubly ionized helium is about 4 M_⊙ below the surface of the star. At this time, the dynamical timescale of the envelope, approximated by

Fig. 2. Restructuring of the envelope during the initial transient. Panel a shows the radius of the envelope (similar to Fig. 1a). The partial ionization zones and the mass coordinates of the recombination fronts of hydrogen and helium are shown in panel b, as well as the total mass of the star. The recombination fronts are defined where the ionization fraction of the respective species reaches 0.5. The partial ionization zones are defined to be the layers where the ionization fraction of the respective species is between 10 and 90%. Panel c shows the specific entropy in units of Avogadro’s number NA and the Boltzmann constant kB, as well as the density. Both quantities are reported at a mass coordinate of 5.5 M⊙, which is close to the bottom of the convective envelope. Panel d shows the radiative cooling timescale in units of the dynamical timescale. The total energy of the convective envelope is shown in panel e. The gray band indicates the catastrophic cooling event at t ≈ 32.5 yr.

$$\begin{array}{c}\hfill {\tau }_{\mathrm{dyn}}=\sqrt{\frac{R^3}{G{M}_{\mathrm{env}}}},\end{array}$$(1)

is longer than the radiative cooling timescale of the envelope (Fig. 2d), approximated by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{KH}}=\frac{GM{M}_{\mathrm{env}}}{2RL}·\end{array}$$(2)

Here, R is the radius of the RSG, M is its mass, G is the gravitational constant, M_env is the mass of the convective envelope, and L is the luminosity of the RSG. Whenever τ_KH < τ_dyn, the envelope is losing energy via radiation faster than it can respond dynamically (c.f. Clayton et al. 2017). This causes a catastrophic cooling event at the surface, where a significant fraction of the internal energy of the envelope is radiated away through the optically thin outer 2 M_⊙ layers. As a consequence, the total energy (sum of potential, internal, and kinetic energy) of the convective envelope decreases (Fig. 2e).

In the outer 2 M_⊙ layers, the opacity is orders of magnitude lower compared to the rest of the envelope, because hydrogen and helium are fully recombined (log(κ/cm² g⁻¹) = 1.8 in the hydrogen partial-ionization zone versus log(κ/cm² g⁻¹) <  − 3 in the neutral layers). Therefore, the neutral layers are stable against convection. As the expanded envelope contracts, the neutral, radiative layers become re-ionized and get mixed with the convective layers further below. This mixes the low-entropy radiative layers with the high-entropy convective layers further inside, lowering the average entropy in the entire convective envelope. Additionally, we observe that the entropy gradient changes during the restructuring event. The entropy gradient becomes less negative, which can be seen by a smaller spacing between the lines of constant mass after the restructuring in Fig. 3. As a consequence of the altered thermal structure of the envelope, the density structure adjusts accordingly. We find on average larger densities deep in the convective envelope and lower densities closer to the surface compared to the original hydrostatic structure. The shallower entropy gradient might also have consequences for the efficiency of the convective energy transport. Such restructuring phases have been previously found and are described in Ya’Ari & Tuchman (1996), Lebzelter & Wood (2005), and Clayton (2018). They seem to be unavoidable for pulsating envelopes once the dynamical and thermal timescales become comparable, even if the initial growth of the pulsations is hindered (see Appendix A).

Fig. 3. Specific entropy, s, in units of Avogadro’s number, NA, and the Boltzmann constant, kB, at various mass coordinates inside the convective envelope. The gray band indicates the restructuring event at t ≈ 32.5 yr.

We expect the catastrophic cooling events would be observable as an infrared transient because of the large amount of energy radiated away in a short duration, and because of the cool effective temperature of the RSGs. For the 15 M_⊙ RSG model, the maximum luminosity during the catastrophic cooling event reaches 7 × 10⁴⁰ erg s⁻¹. In total, 5 × 10⁴⁶ (7 × 10⁴⁶) erg are radiated away during a period of 30 (120) days. Such timescales and energetics are very similar to observed luminous red novae (e.g., Pastorello et al. 2021).

After the restructuring phase, the pulsations continue at smaller but varying amplitude (35 yr ≲ t ≲ 80 yr). During this phase, the envelope is thermally and dynamically adjusting to the structural changes of the catastrophic cooling event. At a time of around 80 yr, the density and entropy converge to a new structure, indicating that the envelope has reached a new thermal and dynamical equilibrium state.

Once the envelope is adjusted to the changes after about 80 yr, the pulsations become periodic. This phase lasts for several hundred years. Because the initial 80 yr of the hydrodynamic simulation is still affected by the switching from the hydrostatic to the hydrodynamic model, it seems unlikely that this phase represents the behavior of the true stellar envelope. In a more realistic scenario, we envision that the envelope starts to pulsate gradually once the star climbs the RSG branch. Thus, it increases its L/M ratio and experiences the envelope restructuring once the amplitudes are large enough for the catastrophic cooling event to set in. The catastrophic cooling event itself is unavoidable because with the growing amplitude of the pulsations, the ratio of the thermal to dynamical timescale decreases, eventually becoming on the order of unity, and triggering a thermal restructuring of the envelope. The exact path of reaching the catastrophic cooling event is expected to occur over a much longer timescale compared to the one we simulate here. Afterward, we find that the steady-state pulsations are identical, even if the evolution and the growth of the pulsations that lead to the catastrophic cooling are changed, for example, by applying a damping force in the envelope (see Appendix A). Hence, we focus on the steady-state pulsations with t > 80 yr in the following analysis of the pulsations.

3.2. Periodic pulsations of RSG envelope

For t > 80 yr, the radial pulsations are periodic with a pulsation period of ${817}_{-11}^{+8}$ d that changes by less than 2% during the simulated time. This period is determined using a Lomb-Scargle periodogram and the uncertainty is the full width at half maximum. For comparison, the initial dynamical timescale is about 230 d. A zoom on one pulsation cycle at 102 yr–104 yr is shown in Figs. 1b and 1c. The luminosity varies by more than one order of magnitude during the pulsation cycle (4.45 ≤ log(L/L_⊙)≤5.90) and the effective temperature varies between 2200 and 6820 K (i.e., log(T_eff/K) = 3.342 and 3.834, respectively). These variations are much larger than typical uncertainties on SN progenitor observations (Δlog(L/L_⊙)≈0.2 and Δlog(T_eff/K)≈0.05; Smartt 2015).

We indicate six points during the pulsation cycle (A–F), which are chosen to highlight the characteristic stages during one pulsation cycle. Point A is chosen to be at the maximum radial extent. Points B and F are chosen to be close to the hydrostatic radius of the RSG at 1024 R_⊙, with point B in the contraction phase and point F in the expansion phase of the pulsation. Points C and E are chosen to be at the local minima of the radial extent of the RSG. These two points are separated by a small expansion phase that lasts less than 100 d. Point D is chosen to be at the local maximum radius between points C and E.

The pulsations are mainly driven by the combination of the κ and γ mechanism (κγ mechanism). The former is connected to the changing opacity in partial ionization layers of hydrogen and helium, and the latter corresponds to the released (consumed) energy by a change in the ionization state of hydrogen and helium during the expansion (contraction). The entire driving mechanism is visualized in Fig. 4, which shows the recombination structure of the RSG at points A–F.

Fig. 4. Recombination structure of the RSG envelope at points A–F. The six colored circles show the radius of the RSG at points A–F. The coloring inside the circles shows the specific released recombination energy ϵrec as logmod(ϵrec/103 erg s−1 g−1) on a linear radial scale, with logmod(x) = sign(x)log10(|x|+1) (John & Draper 1980). Red corresponds to energy released by recombination, while blue corresponds to energy being removed via ionization. The square insets at each point show the recombination luminosity Lrec, i = ∫Menvϵrec, i dm for the species i = H+,  He+ and He2+ on a linear scale. Hatched regions indicate the partial ionization zone for hydrogen and helium, defined where the ionization fraction of the corresponding species is between 1% and 99%. Light gray shading at points A and B corresponds to neutral layers (ionization fractions less than 1%). The central plot shows the radius evolution of one pulsation cycle and indicates the 6 highlighted points A–F. The temporal evolution of the recombination structure is available as an online movie.

At point B, the entire envelope is in a contraction phase and both hydrogen and helium are being re-ionized (blue color in Fig. 4). The process of re-ionization takes energy from the pulsation and converts it to ionization energy (γ mechanism). Additionally, the partial ionization zones of hydrogen and helium increase the opacity, therefore increasing the radiation pressure and Eddington factor, which opposes the contraction of the envelope (κ mechanism; see Appendix E).

As the envelope contracts, hydrogen and helium get completely re-ionized. At point C, the entire envelope is ionized. The outermost layers fall back onto layers of larger density and pressure; they get compressed, bounce off, and re-expand (c.f. Figs. E.1 and E.2). As they cool down during re-expansion, H⁺ and He⁺ start recombining (red color in Fig. 4) and drive the expansion of the surface layers. Only the outer 0.1 M_⊙ are expanding, while the layers further inside the envelope keep on contracting and cause further re-ionization of He²⁺. Because only the immediate surface layers are expanding, the expansion phase starting at point C reaches only a small maximal radial extent at point D of 675 R_⊙. At this point, all the recombination energy from H and He+ of the expanding layers is released. Another contraction phase follows, during which the surface layers are again re-ionized. Meanwhile, layers further inside the envelope keep contracting and re-ionizing He²⁺ until enough pressure (both gas and radiation pressure) builds up to stop the contraction and start the expansion of the interior envelope layers.

At point E, the contracting surface collides with the expanding interior envelope layers, reaching another minimum in total radial extent. This point marks the smallest radial extent throughout the entire pulsation cycle with R = 622 R_⊙. The interior layers are accelerated by the recombination of He²⁺, which marks the main difference to point C, where He²⁺ is being re-ionized.

At point F, the whole envelope expands. The envelope cools down, which causes the recombination of hydrogen and helium. The released recombination energy further accelerates the expansion, driving it to large radii. Most of the released recombination energy is from He²⁺ with smaller contributions of hydrogen and He⁺. The recombination of helium stops already before reaching the maximum extent of the envelope at point A, while the recombination of hydrogen continues for a longer time, giving the surface layers an extra push.

At point A, the star is at its maximal radial extent of 1568 R_⊙. Now, the surface layer starts to contract. However, because of hydrogen recombination, there is a subsurface layer with an outward radial velocity causing a density inversion (Fig. 5). Layers below the hydrogen recombination zone are contracting and heating up, leading to the re-ionization of helium. As a result, the entire envelope contracts, and we are again at the state in point B. In this manner, the pulsation cycle driven by this κγ mechanism is sustained and remains stable.

Fig. 5. Density profiles of the RSG envelope at points A–F during the pulsation cycle in terms of mass coordinate (panel a) and radius coordinate (panel b). The highlighted regions show hydrogen and helium partial ionization zones where the ionization fraction is between 10% and 90%. The density profile of the hydrostatic model is shown for comparison. Arrows along the profiles indicate the radial velocity, both the direction (inwards or negative; outwards or positive) and the magnitude, and highlight layers that are expanding/contracting. The inset in panel a shows the density structure of the interior layers, up to where the core was cut out. The shading indicates the core region of the RSG. The inset in panel b shows the radius evolution of one pulsation cycle and indicates the points A–F. The temporal evolutions are available as an online movie panel a and an online movie panel b.

3.3. Density variations by the pulsations

The large amplitude pulsation of the RSG, driven by the κγ mechanism, significantly alters the density structure of the envelope (Fig. 5). For mass coordinates below 10 M_⊙, we find consistently larger densities compared to the hydrostatic model. This is a consequence of the catastrophic cooling event at t = 32.5 yr that reduces the total energy of the envelope and leads to the rearrangement of the envelope structure. As expected, we find lower densities for the models in the expanded states compared to the models in the compressed states. The largest densities in the deep interior are reached at point D. Although the stellar radius of point D is larger than at points C and E, the deeper layers keep contracting between points C and D, causing the density to increase.

The density near the surface (m > 11.5 M_⊙) shows a lot of structure, and we find large jumps in density. Most of this structure can be directly connected to the recombination and re-ionization processes happening in the envelope. There is a large density inversion at m = 11.7 M_⊙ (r = 1450 R_⊙) in the density profile of point A. Hydrogen recombination at 11.5 M_⊙ accelerates low-density material that crashes into high-density surface layers, which are already contracting. The high-density surface layer is built up by the continuous expansion driven by hydrogen recombination. At point B, there is a large drop in density at 12 M_⊙, which is an immediate consequence of the surface layers falling back onto the deeper, ionized layers. A shock front develops at the interface where the neutral surface layers are re-ionized. The post-shock density is about ten times larger than the pre-shock density (10⁻⁸ g cm⁻³ post-shock versus 10⁻⁹ g cm⁻³ pre-shock). At points C and E, we find smaller density inversions just below the surface. Layers above the density inversion participate in the short-lived expansion phase between points C and E, while layers further below keep contracting and increasing their density. At point F, we find a smooth density profile, as the whole envelope is expanding coherently. However, just below the surface, there is already a pile-up of high-density material from hydrogen recombination. This layer keeps on growing further until point A. Note that the density profiles at points F and B show very different structures, even though the radius of the star is the same.

In the layers below the convective envelope, the density structure is unchanged because the pulsations are confined to the envelope. We find variations of less than 0.5%, which decrease deeper inside the star, indicating that the pulsations are damped in these layers.

3.4. Dust models of RSG pulsations

We used the DM models (see Sect. 2.3) to construct light curves of the pulsating RSG with varying amounts of dust around it. The light curves in the V and K bands are shown in Fig. 6.

Fig. 6. V- and K-band light curves of the pulsating RSG surrounded by a dust shell. The parameter τ0 sets the optical depth at the beginning of the pulsation cycle (ϕ = 0). The light dotted lines indicate regions outside the MARCS grid (too low Teff) that do not allow us to make photometric estimates. For better visualization, two full pulsation cycles are shown. The light, horizontal lines show the time-averaged magnitudes. For comparison, the phase-folded, visual light curve of VX Sgr is shown as well, with data taken between JD = 2640000 and 2644000, and assuming a period of 747 days.

The light curves do not follow a sinusoidal modulation of the brightness, but show a rather complex morphology. We find two short-duration peaks, corresponding to points C and E (i.e., the states of maximum compression). Around the maximum extent of the RSG (point A, ϕ = 0), the light curves vary more smoothly. In the V-band light curves, we find that the brightness is slowly decreasing and might even reach a plateau phase, depending on the amount of dust. In the K-band light curves, this phase corresponds to the maximum light of the star. Unfortunately, we are unable to predict the cool phases of the light curves because the underlying MARCS models are only available for T_eff ≥ 2400 K.

The amount of dust, modeled via different optical depths, τ₀, has a significant influence on the amplitude of the V-band light curve. We varied τ₀ between 0.5 and 10, which corresponds to total dust masses between 10⁻⁴ and 10⁻² M_⊙, in line with expectations of dusty RSGs such as VY CMa (Humphreys & Jones 2022), NML Cyg (De Beck et al. 2025), or the progenitor of SN 2023ixf (Niu et al. 2023). For low amounts of dust, we find a total brightness variation of ≈10 mag, while for increasing amounts of dust, we find a brightness variation of up to 25 mag. The K-band light curves are less influenced by the amount of dust around the RSG. The brightness variation for the K band is around 2.0 to 2.5 mag with increasing amounts of dust.

We find that for the V band, the more dust is present around the RSG, the lower is the brightness at any time during the pulsation. This trend does not hold for the K-band light curve. While the average brightness decreases with increasing amounts of dust, we also find times during the pulsation cycle (e.g., point E at ϕ = 0.55) where the opposite ordering is the case. That is, we find that the RSG is brighter in the model with more dust.

For comparison, we show the visual light curve of VX Sgr⁴, a semi-regular pulsating RSG with a large amplitude in the optical. The peak-to-peak variability of VX Sgr is about 6 − 7 mag, which is lower than the dust free (τ₀ → 0) peak-to-peak variability of about 10 mag of the model. Note that the observations of VX Sgr are visual magnitudes, and the computed light curves are in the V-band. There is also some despute whether VX Sgr is an RSG, an extreme AGB star or even a Thorne-Żytkov object (Tabernero et al. 2021). In any case, the envelope structure should be very similar and comparable to the 15 M_⊙ RSG model. For a more detailed comparison to observations, see Sect. 5.4.

4. Supernova light curves

We modeled the SN light curves for the RSG at various stages throughout the pulsation cycle. Although the pulsation simulations are computed at the end of core carbon burning, we expect the pulsations to continue until the onset of core collapse because the dynamical and thermal timescales of the envelope are larger than the nuclear timescale of the core. Therefore, we do not expect that the timing of the core collapse is in any way related to the pulsation phase of the envelope. This means core collapse can happen at any time during the pulsation cycle of the envelope with equal probability. We highlight the SN light curves corresponding to explosions at points A–F as representative states of the envelope during one pulsation cycle.

4.1. Light curves at different pulsation phases

We simulated the SN explosion of the RSG at different phases of the pulsation with the same explosion parameters, derived by applying the model of Müller et al. (2016) (see Sect. 2.4). The resulting light curves for each pulsation phase are shown in Fig. 7a. We find that the light curve shapes and decline rates are strongly impacted by the pulsation of the RSG. At all phases of the pulsation, the obtained light curves are different from the typical Type II-P shaped light curve obtained by simulating the explosion of a hydrostatic RSG model, which is the standard procedure in this field (see the dash-dotted line in Fig. 7a). At large radii (point A), the light curve is declining fast during the first 80 d (1.1 dex/100 d in L_bol or 3 mag/100 d in the V band⁵), and might have classically been referred to as a Type II-L SN. At later times, the light curve shows the characteristic exponential decay when the radiation output is dominated by the contribution of radioactive ⁵⁶Ni. In contrast, at small radii (points C–E) the light curve has the long (∼80 d) plateau phase, i.e., small decline rates (0.45 dex/100 d in L_bol or 1.37 mag/100 d in the V band), characteristic of II-P SNe, and reaches lower luminosities. At intermediate radii (points B and F), when the star has approximately the hydrostatic radius, the light curve shows an early excess emission followed by a decline with decline rates between the two extremes shown above.

Fig. 7. Bolometric light curves (panel a) and photospheric velocities (panel b) for SNe at points A–F during the pulsation cycle. The explosion of the hydrostatic model is shown for comparison. The inset shows the radius evolution of one pulsation cycle and indicates the points A–F. Once the photospheric temperature drops below 103.75 K, the velocities are no longer reliable because SNEC cannot accurately estimate the location of the photosphere due to limitations in low-temperature opacities (Morozova et al. 2015). The complete ϕ-evolution of the light curve and the photometric velocity is available as an online movie.

We show the corresponding photospheric velocities v_phot in Fig. 7b. There is a clear distinction between points A and F, and all other points. For both points A and F, we find that after the shock passed through the ejecta, v_phot stays constant for up to one month before decreasing. At the other points, v_phot peaks directly after shock breakout and then steadily declines. At these points, the maximum v_phot is also the highest.

Both effects, the higher luminosities and the lower v_phot at points A and F compared to the other points, can be directly connected to the different envelope structures at the time of explosion. At points A and F, the envelope is extended to beyond 1000 R_⊙. This leads to more deceleration compared to the compact progenitors, i.e., lower photospheric velocities. The higher the deceleration, the more explosion energy is converted to radiation, which leads to a higher luminosity. At point B, the radius of the progenitor is the same as point F, but the density structure is different, having larger densities in the outer layers of the envelope (Fig. 5). This results in a different light curve and photospheric velocity at point B compared to point F. In summary, we find that taking the hydrodynamical changes of RSGs into account, i.e., their radial pulsations driven by the κγ mechanism, naturally reproduces some of the diversity of hydrogen-rich SN light curves (Anderson et al. 2014).

4.2. Connecting light curve features to the progenitor structure

As was noted in the previous section, the light curves we obtained at different pulsation phases have different shapes. They also display features, such as the short plateau-like feature of about 10 d at point A. We investigate the origin of these shapes and features by studying the connection to the progenitor structure in Fig. 8; in panel b of this figure, we present the pre-SN density profiles at different phases of the evolution. These are highly affected by the ionization and recombination states described in Sect. 3.3. As expected, at the point of the largest radius expansion (point A), the stellar structure has a particularly low density in the envelope, while at the point when the smallest radius is reached, it has the highest overall density. The density profiles show prominent density inversions near the surface that are directly related to the ionization structure (see Sect. 3.3).

Fig. 8. Connecting feature in the SN light curves to features in the density profiles of the RSG envelope. Panel a is similar to Fig. 7a and panel b is similar to Fig. 5. The density profiles are shown as a function of the mass until the surface. Colored markers connect features in the density profiles to features in the SN light curves. Features in the light curve arise when the SN photosphere reaches the mass coordinates highlighted in the density profile.

While the SN shock goes through the stellar structure, it heats and accelerates the layers it encounters. Once it has traversed most of the entire stellar structure and the first photons escape, the SN light curve begins. For different pulsation phases, the progenitor size varies greatly. This is reflected by the time it takes until the peak luminosity is reached, which increases for larger RSG radii (light curves A vs. C in Fig. 8).

The SN shock leaves the ejecta completely ionized, and the opacity increases greatly in the layers that were originally recombined. Photons are trapped behind the optically thick layers. As a cooling mechanism, the ejecta expands rapidly and the temperature and luminosity drop, leading to a drop in the light curves. At this point, the photosphere remains at the mass coordinate of the surface of the pre-explosion star (c.f. Fig. F.1). The rate at which the luminosity drops reflects the state of the progenitor. The more expanded the outer layers, the lower the density, and the slower the initial decline rate of the light curve.

Once the hydrogen-rich ejecta have expanded and cooled sufficiently for the recombination of hydrogen atoms to occur, the photosphere moves inward, and the photons trapped behind it are released, contributing to the SN light curve. In addition, the temperature decreases less rapidly. By following the location of the photosphere during the explosion, the features in the SN light curves can be traced back to the moments when the photosphere reaches certain parts of the stellar layers. We highlight these moments by markers on the density profiles and SN light curves in Fig. 8.

We find that when the photosphere reaches the layers where the density inversions connected to the partial-ionization zones occur, features in the light curves arise. At point A, there is a density inversion in the pre-SN density profile at log(M_⋆ − m)/M_⊙ = −0.3. After the SN, the photosphere reaches this particular mass coordinate at a time of 30 d. At this particular time, we see a short plateau feature in the corresponding light curve (Fig. 8). The density inversion at point A arises from the released recombination energy of hydrogen (see Sect. 3.3). Similarly, at point B we find that the density increases rapidly by one order of magnitude at log(M_⋆ − m)/M_⊙ = −0.6 because of the re-ionization of hydrogen and helium. This particular mass coordinate is reached by the photosphere 25 d after the explosion and marks the transition of the light curve from the near-linear decline to the plateau phase. A similar feature is visible in the light curves corresponding to point F, where the decline rate changes about 20 d after the explosion. This correlates with the density inversion found at log(M_⋆ − m)/M_⊙ = −1.1, caused by hydrogen recombination. Point D behaves qualitatively similar to point F. The density inversion at log(M_⋆ − m)/M_⊙ = −1.5 is caused by hydrogen recombination during the intermediate expansion phase. When the photosphere reaches this mass coordinate at about 18 d after the explosion, the decline rate changes (i.e., from linear decline to plateau). Points C and E show only weak density inversions that get washed out by the SN shock. When the photosphere passes these mass coordinates, no significant changes in the light curves appear.

5. Discussion

5.1. Timing of RSG pulsations

It has been shown by Heger et al. (1997), Yoon & Cantiello (2010), and also Clayton (2018) that the growth rate of the pulsations increases with an increasing luminosity-to-mass ratio, L/M, which is closely connected to the Eddington factor. Therefore, we expect stronger pulsations as stars evolve in the RSG phase because they typically get brighter while losing mass via winds. For the 15 M_⊙ RSG model that we analyzed, we did not find any pulsations during core helium burning. The pulsations only become significant once the star evolves to burn carbon in its center. Based on the work by Heger et al. (1997) and Yoon & Cantiello (2010), weaker pulsations are expected in lower mass RSGs to occur only shortly before core collapse. For higher mass RSGs, the pulsations are expected to start earlier in the evolution and become more extreme.

Clayton et al. (2017) showed that similar pulsations might also be present in the context of common-envelope evolution. They find that depending on the energetics of the common-envelope event, the pulsations might even lead to the dynamical ejection of the outer envelope.

The pulsation mechanism itself depends only weakly on the metallicity of the star, because the partial ionization zones of hydrogen and helium are mainly responsible for driving the pulsation. Therefore, RSGs of similar mass as considered here but with a different metallicity should behave similarly regarding the pulsations. However, the mass loss rate during earlier evolutionary phases might heavily depend on the metallicity (e.g., Puls et al. 2008; Mokiem et al. 2007). In such cases, the mass loss has typically only a minor effect on the luminosity of the star. Hence, the L/M ratio of higher metallicity stars is larger, and we expect even stronger pulsations.

For the 15 M_⊙ RSG, a SN with a fast-declining light curve (classically Type II-L) is more likely compared to a slow-declining (Type II-P) light curve, because the RSG spends more time in the extended phase. Observationally, many more II-P SNe are found compared to fast-declining (II-L) SNe amongst all core-collapse SNe – Smartt (2009) report 59% II-P and 3% II-L while Smith et al. (2011) report 48.2% II-P and 6.4% II-L. Therefore, it might seem as if our model over-predicts the number of fast-declining SNe. However, because of the initial-mass function (e.g. Salpeter 1955; Kroupa 2001), lower-mass RSGs are more common SN progenitors compared to higher-mass RSGs. As described above, we do not expect lower-mass RSGs to pulsate significantly at core collapse Paper II. Hence, we would expect most lower-mass RSGs to explode as Type II-P SNe only. Extending on this, we predict that the progenitors of fast-declining (Type II-L) SNe that originate from this mechanism are all massive RSGs.

5.2. Mass loss during the RSG phase

Inferring the winds of RSGs either empirically or via models is challenging. This has led to a plethora of RSG wind prescriptions (de Jager et al. 1988; Nieuwenhuijzen & de Jager 1990; van Loon et al. 2005; Beasor et al. 2020; Kee et al. 2021; Antoniadis et al. 2024; for a review, see Smith 2014). These prescriptions can differ by more than one order of magnitude. In our models, we assume that the RSG wind mass loss follows the prescriptions of de Jager et al. (1988). Varying the wind prescription in the RSG phase might have dramatic effects on the pulsations of the envelope. For higher (lower) winds, the total convective envelope mass before the onset of core collapse can be lower (higher). This can impact the density structure of the envelope and L/M, which in the end are responsible for the exact morphology of the pulsation and possibly also the pulsation period (Gough et al. 1965; Heger et al. 1997).

It has also been suggested that the pulsations themselves are responsible for the high mass loss rates of RSGs. Yoon & Cantiello (2010) suggest that the pulsations can drive a superwind during the final evolution as an RSG. These winds would increase as the pulsations become stronger, i.e., in more massive and more evolved RSGs. Alternatively, dynamical mass ejections remove some envelope material on a short timescale (Tuchman et al. 1978; Clayton et al. 2017; Clayton 2018). Observationally, the highest mass loss rates in RSGs are found for objects like VY CMa (Humphreys & Jones 2022), VX Sgr (Decin et al. 2024), and NML Cyg (De Beck et al. 2025), which all show significant pulsations (e.g., Monnier et al. 1997; Kiss et al. 2006).

In either case, the enhanced mass loss leads to a larger L/M ratio, because the luminosity at the RSG stage is only determined by the carbon-oxygen core mass (e.g., Temaj et al. 2024). This leads to more extreme pulsations that continue to remove the envelope. Depending on the overall strength of such winds or mass loss, the RSG might be able to evolve again toward hotter temperatures once the mass of the hydrogen-rich envelope decreases below some critical value. Then, the envelope is no longer unstable to pulsations, and the star might end its life as a yellow supergiant. In this case, an explosion as a Type IIb or even Type Ib is expected.

For simulations of the 15 M_⊙ RSG presented above, we turned off any mass loss during the dynamical evolution. Because the luminosity and radius of the star change by up to one order of magnitude during one pulsation cycle, any wind determined by classical wind mass-loss prescriptions (e.g., de Jager et al. 1988) would be modulated directly by the changes in these surface properties. However, it is unclear whether such prescriptions provide a meaningful wind mass-loss for dynamical models. Therefore, we cannot make statements about pulsationally driven superwinds, but refer to the results obtained by Yoon & Cantiello (2010). We do not find any dynamical mass ejection for the 15 M_⊙ simulations, but expect these to set in for more massive RSG with larger L/M. The overall effect of successive mass ejections might then be similar to a pulsationally driven superwind, producing a shell-like structure of CSM around the RSG, similar to the CSM observed around VY CMa (Singh et al. 2023) and NML Cyg (Andrews et al. 2022), or thermally pulsing asymptotic giant branch stars (e.g., Olofsson et al. 1988; Mattsson et al. 2007). Additionally, such mass loss should only weakly depend on the metallicity of the star, because the pulsation mechanism described in Sect. 3.2 does not rely on the metal content in the envelope of the RSG.

In any case, all these effects might leave some CSM around the SN progenitor. This is true especially for more massive RSGs, as they exhibit the strongest pulsation, in which we expect dynamical mass ejections (Clayton 2018). This has direct implications for the missing RSG problem (Smartt et al. 2009; Smartt 2015) and could explain the maximum luminosity observed in RSGs (Davies et al. 2018).

The SN ejecta are also expected to interact with the CSM that is produced by enhanced mass loss during the pulsation. The CSM interaction typically leads to excess emission during the early light curve (Moriya et al. 2011), as observed in many Type II SNe (Anderson et al. 2024). Depending on the exact nature of the CSM, the star might even produce Type IIn SNe, i.e., hydrogen-rich SNe with narrow line emission features in the spectra. Including the CSM effects on the light curves is beyond the scope of this paper, but is discussed in detail in Paper II.

5.3. Pulsating RSGs in binary stars

If a pulsating RSG is part of a binary star system, many evolutionary scenarios can open up, depending on the initial separation (e.g., Podsiadlowski et al. 1992, 1993; Nomoto et al. 1993, 1995; Ercolino et al. 2024; Dessart et al. 2024). There could be episodic mass transfer to the companion star, possibly leading to a complex and shell-like CSM structure (Landri & Pejcha 2024). The evolution of the binary orbit is most likely very sensitive to the initial conditions of such an episodic mass transfer phase, and it is plausible that such a phase might increase or decrease the eccentricity of the binary depending on how the pulsation period relates to the orbital period and how the pulsation period changes upon mass transfer. Naively, one expects the pulsations and orbital timescale to be similar upon mass transfer, as the pulsations are a dynamical phenomenon and the dynamical timescale is similar to the orbital period⁶. However, the interaction between mass transfer and the pulsations may result in a nontrivial evolution of the system.

5.4. Comparison to observed pulsating RSGs

A large number of observed RSGs show period variability (Kiss et al. 2006). Based on theoretical arguments, the pulsation period of RSGs is supposed to scale linearly with L/M (Gough et al. 1965). From observations, period-luminosity (PL) relations can be measured to use RSGs as standard candles for distance estimation (e.g., Jurcevic et al. 2000). Practically, the absolute magnitude in the K band, M_K, is used because it is more accessible than the bolometric luminosity, and it is only little affected by dust attenuation. Based on the models of the 15 M_⊙ RSG with no dust, we find a period of ${817}_{-11}^{+8}$ d and ⟨M_K⟩=−10 mag (Fig. 6b). Chatys et al. (2019) used RSGs in the Galaxy and the Large Magellanic Cloud to derive a PL relation. For a period of 817 d, the PL relation would expect M_K to be between −9.5 mag and −11.5 mag. Using the PL relation from Soraisam et al. (2018) based on RSGs in M31, we would expect M_K ≈ −11.2 mag with a spread of about 0.5 mag. Lastly, using the PL relations from Ren et al. (2019), derived from RSGs in M33 and M31, one finds M_K = −11.1 ± 1.4 mag (M33) and M_K = −11.1 ± 1.1 mag (M31). These expectations from PL relations agree within the quoted uncertainties with our models.

Observed amplitudes of the K-band variability reach up to 0.5 mag and tend to increase with the pulsation period (Yang et al. 2018; Chatys et al. 2019; Soraisam et al. 2018; Ren et al. 2019). This trend is also suggested by theoretical models (Heger et al. 1997; Yoon & Cantiello 2010). However, we do find larger amplitudes in our models, independent of the amount of dust around the RSG. The amplitudes we find in the V band for models without any dust are ΔM_V ≈ 5 mag and reach values of > 15 mag for considerable amounts of dust (high τ₀). Some observed RSGs reach ΔM_V ≈ 2 − 3 mag (e.g., S Per, VY CMa; Kiss et al. 2006), with VX Sgr showing the largest brightness variations of up to 7 mag. More typically, amplitudes of 1 mag are found in the V band (Levesque et al. 2007). There are multiple reasons for the discrepancy of the amplitude from our models with observation, such as the treatment of convection in 1D models (see Sect. 5.6). Additionally, our 1D models only allow for radial pulsations, but non-radial pulsations might be present in the envelopes of RSGs. These could lead to a decrease in the amplitude of the radial pulsations. Most importantly, the observed RSGs are all expected to be core-helium burning because of the much shorter timescale of any later burning phases. We only find these pulsations in the very late phases of evolution, i.e., core carbon burning and later. We expect that most of the observed RSGs in the Galaxy or nearby galaxies are currently in the core-helium burning phase because of the much longer timescale. The post-core carbon burning phase until core collapse lasts typically about 10² yr, which is about 0.01% of the entire RSG phase. Moreover, assuming that the pulsations only happen in the more massive RSGs (> 12 M_⊙), there is a chance of less than 1 : 10⁴ of catching an RSG in the phase in which it is pulsating because of an unstable envelope. Because only a couple of 10³ RSGs have been monitored closely in the Milky Way and nearby galaxies (Dai et al. 2025a,b), we expect that zero to one RSGs could have been caught with such large amplitude pulsations. Therefore, we do expect larger amplitudes in our models compared to the observations (see Sect. 5.1). We predict that if observed, such rare large-amplitude pulsations would be an indication of the imminent collapse of such stars. Massive RSGs that already pulsate with a large amplitude during the core-helium burning have likely already lost a significant amount of their envelope via dynamical mass ejections (Clayton 2018) or an enhanced wind (Yoon & Cantiello 2010) and no longer appear as an RSG, linking to the missing RSG problem.

Pre-explosion imaging of the progenitors to SN 2023ixf and SN 2024ggi show signs of periodical variability (e.g., Jencson et al. 2023; Kilpatrick et al. 2023; Soraisam et al. 2023; Xiang et al. 2024a,b). In the case of SN 2023ixf, a dense and close CSM can be inferred from early flash spectroscopy (e.g., Li et al. 2024; Zimmerman et al. 2024). Discussing the details of SN 2023ixf and SN 2024ggi within the context of our model of pulsations RSG and their impact on light curves is beyond the scope of this paper, but will be analyzed and discussed in detail in Paper II.

5.5. Implications of envelope restructuring

The envelope restructuring, observed at around t = 32.5 yr, immediately changes the density structure of the envelope and deviates from the original, hydrostatic envelope structure by up to one order of magnitude. This also means that the pulsations before and after the restructuring event can change significantly, including the pulsation period (Ya’Ari & Tuchman 1996; Tuchman et al. 1978). Predicting post-restructuring pulsation properties from a hydrostatic structure before the restructuring based on linear analysis is impossible.

Usually, the hydrostatic stellar structure at the onset of core collapse is used in simulating SN light curves. However, if such an envelope restructuring occurs before the onset of core collapse, hydrostatic models are unable to capture this effect, independent of how large the amplitude of the pulsations is.

Goldberg et al. (2020) model the fundamental and first overtone pulsation of an initially 18 M_⊙ RSG just before the onset of core collapse, and then studied their effect on the resulting SN light curve. They find that the fundamental-mode pulsations do not significantly modify the SN light curve in unusual ways. The pulsations decrease (increase) the density in the entire envelope upon expansion (compression). These changes in the SN light curve from the fundamental-mode pulsations can be recovered with typical scaling relations (e.g., Popov 1993; Goldberg et al. 2019) by considering a progenitor with the same mass but larger (smaller) radius. They do report that the amplitude of the pulsations is still growing at the onset of core collapse. However, they do not find any restructuring of the envelope. This is most likely because they start simulating the pulsations at most 1750 d (just over 3 full pulsation cycles) before the core collapses. This means that the pulsations did not have enough time to grow in amplitude such that a catastrophic cooling event can be triggered, which ultimately causes the envelope to restructure. The amplitudes of their pulsations are motivated by the observed brightness variations of RSGs. However, we stress that most observed pulsating RSGs are not beyond core-carbon burning (see Sect. 5.1), and, therefore, we do expect larger amplitudes at core collapse.

In contrast to the fundamental-mode pulsations, Goldberg et al. (2020) find a significant deviation of the SN light curve compared to the hydrostatic progenitor when considering first-overtone pulsations. In this case, the density in the inner part of the envelope increases (decreases) while the density in the outer part of the envelope decreases (increases) upon expansion (compression), with one nodal point in between. In the resulting SN light curve, Goldberg et al. (2020) find an increased brightness at earlier times for an expanded progenitor. This is very similar to our findings. However, in our case, the density in the inner part of the envelope increases because of the restructuring after the cooling catastrophe. The resulting effect on the SN light curves is very similar because the density structure is modified similarly. In our case, we find more severe changes in the SN light curves because the envelope restructuring together with the large amplitude pulsations introduce a larger deviation from the hydrostatic structure compared to the effect of the first-overtone pulsations reported in Goldberg et al. (2020).

In general, simulations of the core evolution and the dynamical envelope evolution at the same time to capture the onset of the pulsations and to simulate them to core collapse are not possible. The difference in the dynamical timescale of the envelope and the nuclear timescale in the core is too large. Nonetheless, it is important to start simulating the pulsations early enough before the core collapses, such that they can fully develop and possibly reach an equilibrium state, because we find that their effect on the final stellar structure and explosion properties is significant.

5.6. Convection in pulsating envelope

The envelopes of RSGs are unstable against convection (i.e., the bulk motion of matter). The direct interplay between convection and pulsation is not well understood (e.g., Wood 2007; Houdek & Dupret 2015; Freytag et al. 2017; Goldberg et al. 2022a; Ahmad et al. 2023; Ma et al. 2024). Convection itself is a multidimensional phenomenon that might trigger non-radial pulsations, especially given that only a few convective cells are expected to be present in the envelope of RSGs (Goldberg et al. 2022a; Ahmad et al. 2025). These can remove energy from the radial pulsations and reduce their amplitudes, weakening the overall effects of the pulsations and potentially smoothening the RSG light curve (Fig. 6). Olivier & Wood (2005, 2006) include turbulent viscosity when modeling pulsations of red giant stars. This introduces damping of the pulsations and reduces their amplitudes. In our models, we find that increasing the mixing-length parameter of convection has a similar effect (see Appendix C). However, the physical mechanisms at work are different because the increased mixing length leads to more efficient energy transport and more compact RSGs. There also exist alternative treatments to classical mixing-length theory (Prandtl 1925; Cox & Giuli 1968), such as time-dependent models for convection (e.g., Kuhfuss 1986; Xiong et al. 1997), which might offer a better description of convection in stellar envelopes. Including such models in simulations of oscillating red giants, Xiong et al. (1997) concluded that turbulent pressure can have an important effect on the pulsations. However, it has been shown by different studies that modifying the treatment of convection has negligible effects on the pulsations; that is, by introducing a phase shift between the pulsations and the convective flux (Langer 1971) or completely freezing the convective flux (Heger et al. 1997). In general, 3D simulations of convective envelopes (e.g., Goldberg et al. 2022a; Freytag & Höfner 2023; Ma et al. 2024) can help us to better constrain the interplay between convection and pulsations.

Recent observations of asymptotic-giant-branch stars have revealed that the size of the convective cells changes throughout the pulsation cycle (Rosales-Guzmán et al. 2024). This might open up a new pathway to directly study the connection between convection and pulsations without the need for numerical simulations.

The treatment and efficiency of convection can directly impact the amplitude and period of the pulsations. This study aims to show that the pre-SN pulsations of RSGs exist and that they have a direct impact on the SN light curve. The role of convection can weaken the results on the SN light curves via lower amplitude pulsations. Nonetheless, we expect the amplitude to increase with the mass of the RSG; therefore, we still expect to see similar effects in higher-mass RSGs.

5.7. Limitation of pulsations models

There are several limitations to the models presented above. They mostly originate from the general limitation of 1D models compared to 3D models.

The envelopes of RSGs are dominated by large convective cells (Goldberg et al. 2022a), which cannot be treated in 1D (see discussion above; Sect. 5.6). By assuming spherical symmetry for 1D simulations, any inhomogeneities in the envelope that can be induced by the large-scale convective motion in 3D are removed. Full 3D radiation-hydrodynamical simulations (e.g., Goldberg et al. 2022a; Ahmad et al. 2023) show that these variations can be significant, especially in the outermost layers of the stars. Consequently, the energy transport by radiation through these layers is changed and can be more efficient than expected by 1D models with standard MLT treatment Goldberg et al. (2022a). This might further have implications for the expected wind mass loss, which could be enhanced in certain locations of the surface, similarly, for example, to the great dimming of Betelgeuse (Dupree et al. 2022). All of these effects can decrease the strength of the pulsations and the resulting diversity in SN properties. The non-sphericity can immediately impact the early evolution of the SNe light curve; the shock breakout will happen at different radii and at different times, lowering the overall luminosity during this phase (Goldberg et al. 2022b).

The structure of the 1D envelope of the RSG shows several prominent density inversions (see Fig. 5). In a 3D realization of such an envelope, one might suspect that these inversions become weaker because they are Rayleigh-Taylor unstable. However, the dynamical timescale is similar to the evolutionary timescales of these features, making it unclear whether there is enough time to completely dissolve them in the envelope structure. Additionally, these density inversions are direct consequences of the recombination processes in the envelope. The released recombination energy during the expansion can continue to drive and maintain such density inversions during the period of the recombination. Any of these effects can result in the reduction of the pulsation amplitudes and reduce the diversity found in the SN properties. These might also smooth the progenitor light curve, removing any short-duration peaks and producing a more sinusoidal light curve (see Fig. 6).

Another limitation is that we need to assume a mixing length parameter α_MLT for the treatment of convection. For our fiducial model, we assume α_MLT = 1.8, which is inspired by calibrations to the Sun but also reproduces the HRD locations of Galactic RSGs (e.g., Ekström et al. 2012). A comparison between observed SN properties and the predictions from RSG progenitor models led to the suggestion of using a larger mixing length parameter (α_MLT ≈ 3) that produces more compact RSG progenitors (e.g., Dessart et al. 2013; Paxton et al. 2018). Three-dimensional radiation-hydrodynamical simulations of convection in RSGs also support α_MLT ≈ 3 − 4 deep in the envelope, but find lower α_MLT near the surface (Goldberg et al. 2022a; Ma et al., in prep.). The usage of MESA’s MLT++ in our simulation reduces the superadiabacity in radiation-dominated convective layers (Paxton et al. 2013), making convection more efficient, as is suggested by the 3D simulations.

Progenitors with more efficient convection (e.g., α_MLT = 3) show at a fixed mass lower amplitude pulsations and less diversity in the SN light curves (see Appendix C). However, when adopting a larger α_MLT, more massive RSGs are still expected to exhibit larger amplitude pulsations. Although the radii of RSGs with a higher α_MLT are generally smaller, they increase with increasing stellar mass, meaning that also the pulsation amplitudes increase. A more massive RSG with a higher α_MLT is expected to show pulsations of similar amplitude as a lower mass RSG with a lower α_MLT. Therefore, increasing the mixing length parameter has the effect of increasing the initial mass of stars at which we expect large-amplitude pulsations to be present.

5.8. Inferring explosion and pre-supernova properties

There exist numerous techniques to estimate explosion and progenitor properties of Type II SNe from the light curve properties (e.g., Popov 1993; Goldberg et al. 2019). We show in Sect. 4 that the explosion of the same star with identical explosion properties but taken at different pulsation phases can lead to significantly different SN light curves. Because the density structure is nontrivially modified by the pulsations compared to the hydrostatic structure, scaling relations fail to connect the SN progenitor properties to the SN light curve (Goldberg et al. 2020). Therefore, inferences of SN progenitor properties from the SN light curve need to be treated with care when the progenitor star is pulsating, especially when the amplitude of the pulsation is large.

Moreover, inferring progenitor properties based on archival imaging can lead to incorrect luminosity and temperature estimates, especially when only one epoch of observations is available. The uncertainties we expect from the trace of the RSG in the HRD during the pulsations can be even larger than observational uncertainties (see Fig. 1c). Therefore, observing the progenitor at just one instance of time will lead to a significant over- or underestimate of its luminosity and/or temperature. Future surveys, such as the Legacy Survey of Space and Time (LSST) at the Vera C. Rubin Observatory (Ivezić et al. 2019), can increase the number of progenitors and may provide even full progenitor light curves.

6. Conclusion

In this paper, we study the pulsations of RSGs and their immediate implications on the light curves of Type II SNe. We follow the final evolutionary stages of a 15 M_⊙ RSG. Our results can be summarized as follows:

• By direct hydrodynamical simulations of the RSG, we find radial pulsations during core carbon burning. During core helium burning, the RSG envelope remains stable. These findings confirm earlier work by Heger et al. (1997), Yoon & Cantiello (2010), and Clayton (2018).

• We show that these pulsations are driven by the κγ mechanism through the cyclic recombination and ionization of hydrogen and singly and doubly ionized helium.

• The pulsation period agrees with the PL relation of observed RSG.S The amplitude of the brightness variations is larger than in observations. We predict that such large-amplitude variations are rare because they only occur shortly before core collapse.

• We find that taking the radial pulsations of massive RSGs into account leads to light curves that are significantly different from light curves obtained for a hydrostatic model.

• The light curves at different pulsation phases naturally reproduce some of the diversity of hydrogen-rich SN light curves. We find different decline rates in the SN light curve, with the fastest-declining light curves for the most extended progenitors and slow-declining light curves from the more compact progenitors. The light curves also show early excess emission without the need to invoke extra CSM around the SN.

• There is a clear distinction between the photospheric velocity evolution of the fast- and slow-declining SNe. For the fast-declining (Type II-L) SNe, the photospheric velocities plateau for up to 30 d and reach the smallest maximum velocities. Slow-declining (Type II-P) SNe reach the larger maximum photospheric velocities that monotonously decrease afterward.

• We find that the light curve shapes and features can be connected to the pre-SN density profiles, which are directly impacted by the locations of the partial-ionization zones.

• The radial pulsations before the SN explosion cause variations in the luminosity and effective temperature of the RSG by up to an order of magnitude. Hence, the locations of massive RSGs in the HRD undergo significant changes before core collapse, implying that only limited constraints can be gained from the observation of pre-SN luminosities of RSGs. This uncertainty is expected to increase for more massive progenitors, where we even expect dynamical mass ejection (Clayton 2018), directly impacting the missing RSG problem.

We aim to extend this study by analyzing a larger parameter space of RGSs; for example, by studying RSGs with different initial masses. This way, we can learn about the pulsations with different amplitudes and their effect on the SN light curves Paper II.

Data availability

Movies associated to Figs. 4, 5, 7 are available at https://www.aanda.org

Acknowledgments

We thank J. Spyromilio for interesting discussions at the MIAPbP 2023 (Interacting Supernovae) that triggered this work. We thank the referee for useful and constructive feedback that helped to improve the paper. We thank J. Goldberg and E. Farag for constructive discussions that helped improve the quality of this work. VAB, EL, FRNS, and PhP acknowledge support from the Klaus Tschira Foundation. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 945806), and is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1-390900948 (the Heidelberg STRUCTURES Excellence Cluster). VAB acknowledges support from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). EL acknowledges support through a start-up grant from the Internal Funds KU Leuven (STG/24/073) and through a Veni grant (VI.Veni.232.205) from the Netherlands Organization for Scientific Research (NWO). We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. Software: astropy (Astropy Collaboration 2013, 2018, 2022), CMasher (van der Velden 2020), matplotlib (Hunter 2007), numpy, (Harris et al. 2020), PyAstronomy (Czesla et al. 2019), PYPHOT (Fouesneau 2025), scipy (Virtanen et al. 2020), TULIPS (Laplace 2022).

References

Appendix A: Damping the pulsations

When initializing the hydrodynamical simulation, spurious velocities might be introduced that can dominate the hydrodynamical evolution of the model. To test whether this is the case, we apply an initial damping force that is slowly reduced over time. This is achieved by extending the momentum equation via

$$\begin{array}{c}\hfill \dot{v}=\gamma v,\end{array}$$(A.1)

with a damping factor γ. This new term causes an exponential decay of any velocities on a timescale of 1/γ. We follow the description of Pakmor et al. (2012) and Ohlmann (2016) and choose an initial damping timescale of τ_dyn/10 that is then exponentially increased to τ_dyn, prescribed by

$$\begin{array}{c}\hfill \gamma (t)=\{\begin{array}{cc}\frac{10}{{\tau }_{\mathrm{dyn}}}& t<10{\tau }_{\mathrm{dyn}}\hfill \\ & \\ \frac{10}{{\tau }_{\mathrm{dyn}}\times {10}^{\frac{t-10{\tau }_{\mathrm{dyn}}}{10{\tau }_{\mathrm{dyn}}}}}& 10{\tau }_{\mathrm{dyn}}\le t\le 20{\tau }_{\mathrm{dyn}}.\hfill \\ & \\ 0& t>20{\tau }_{\mathrm{dyn}}\hfill \end{array}\end{array}$$(A.2)

This ensures that any pulsations caused by the changing damping factor have enough time to decay. After 20τ_dyn, the damping is turned off and the model is allowed to evolve freely.

When applying this damping scheme to the 15 M_⊙ RSG model, we find that the pulsations are quenched for about 10 yr or about 20τ_dyn (Fig. A.1). Then, the model starts to pulsate and experiences a catastrophic cooling event that is accompanied by the restructuring of the envelope. Afterward, we find steady pulsations with an identical pulsation period of ${817}_{-10}^{+8}$ d and the radius evolution differs by less than 0.4%.

Fig. A.1. Comparison of the growth of the pulsations (panel a) and the steady state pulsations (panel b) between the undamped and damped model. The shading during the initial ∼10 yr indicates the magnitude of the damping factor γ.

Once the damping is turned off after about 10 yr, the envelope starts to contract immediately. Afterward, the pulsations start at a much larger amplitude compared to the undamped case, and the catastrophic cooling is already reached after about 24 yr (Fig. A.2). The reason for the more violent initial evolution of the model with damping is that the thermal timescale of the envelope (∼11 yr) is similar to the damping duration. This means that any thermal timescale evolution of the envelope, for example caused by the net energy loss at the surface via radiation, cannot be compensated for by contraction or expansion (Fig. A.2b). Only once the damping term is small enough or turned off entirely does the envelope adjust. In our case, this means that the envelope contracts on the dynamical timescales and initiates the large-amplitude initial pulsation. Nonetheless, we find a growth rate of η ≈ 1.9, which is very close to the undamped case. This indicates that even though the start of the pulsations is more violent, the growth is identical to the undamped case.

Fig. A.2. Restructuring of the envelope for the model with initial damping, similar to Fig. 2. Panel (a) shows the specific entropy in units of Avogadro’s number NA and the Boltzmann constant kB, as well as the density. Both quantities are reported at a mass coordinate of 5.5 M⊙, which is close to the bottom of the convective envelope. The total energy of the convective envelope is shown in panel (b). For comparison, the undamped model is shown as well (same as Fig. 2c and e).

During the catastrophic cooling event, we find that more energy of the envelope is radiated away compared to the undamped case. This results in a different envelope structure immediately after the catastrophic cooling. However, at t ≈ 70 yr, both models converge to the same structure on the thermal timescale (Fig. A.2). This shows that even though the evolution and the growth of the pulsations between the two models differ, they reach the same envelope structure after the catastrophic cooling and the accompanied restructuring.

Generally, the damping has negligible effects on the steady-state pulsations that we are interested in and which can be present for a sustained duration before the SN. This shows that these pulsations are robust, even if the initial evolution and growth of the pulsations are altered.

Appendix B: Different core masses

For our fiducial model of the 15 M_⊙ RSG, we removed the core at a mass coordinate of 4.6 M_⊙. To test whether this choice has a significant impact on the steady-state pulsations, we computed models with a removed core of 4.4, 4.5, and 4.7 M_⊙. The growth of the pulsations up to the catastrophic cooling event and the envelope restructuring are identical between all the models (Fig. B.1). We find differences afterward when the models adjust to the new structure, and the pulsations are not yet periodic. However, once the steady-state pulsations are reached, we find that the pulsation periods differ by less than 1% and the radius evolution differs by less than 0.5%. Therefore, we conclude that the exact location of the core removal has only a minor effect on the steady-state pulsations of the RSG envelope.

Fig. B.1. Comparison of the growth of the pulsations (panel a) and the steady state pulsations (panel b) between the model with different cut-out core masses. The total helium core mass is 5.1 M⊙, defined by a hydrogen abundance less than 0.1.

Appendix C: Dependence of the pulsations on the convective mixing length parameter

In Fig. C.1, we show a comparison of the pulsations of the 15 M_⊙ RSG model (i.e., mixing length parameter α_MLT = 1.8) and a 15 M_⊙ RSG with α_MLT = 3.0. For this, we take the initial stellar structure at the end of core carbon burning and adjust α_MLT to 3.0 in the envelope only. In particular, this means that the core structure is unchanged and L/M is the same, ensuring a meaningful comparison between the models. The larger α_MLT, the more efficient the energy transport via convection. The computational methods are identical to the description in Sects. 2.1 and 2.2. The radius of the RSG with the increased value of α_MLT is 740 R_⊙ compared to 1024 R_⊙ with the lower value of α_MLT. Moreover, we find both a shorter pulsation period (545 ± 5 d) and a lower pulsation amplitude for the high α_MLT RSG. This shows that the treatment of convection directly modulates the strength of the pulsations.

Fig. C.1. Comparison of the pulsations of two RSG models with different αMLT values over two pulsation cycles.

Fig. C.2. Similar to Fig. 7a, but for the progenitor with αMLT = 3, where different colors correspond to different explosion phases. The light curves for the progenitor with αMLT = 1.8 are shown with gray colors as a comparison.

Additionally, the model with the larger α_MLT avoids a catastrophic cooling event, because it stays generally more compact, which causes the thermal timescale (∝R⁻¹) to be always much larger than the dynamical timescale (∝R^3/2). Therefore, no significant envelope restructuring occurs in this model.

We computed SN light curves at characteristic points during the pulsation cycle as described in Sect. 2.4. The explosion parameters are unchanged because the core structure is identical to the α_MLT = 1.8 model. The light curves of the resulting SNe are shown in Fig. C.2. The diversity of the light curves because of the pulsating progenitor is weaker compared to the α_MLT = 1.8 case. This is expected because the amplitude of the pulsation is lower, which causes weaker structural changes in the envelope.

We note that the choice of α_MLT in the convective envelope of RSGs directly impacts the mass loss because of the different radii. This can then lead to different envelope masses at core collapse and therefore also to different L/M. Because L/M directly impacts the pulsation properties, we expect larger differences between the pulsations when α_MLT = 3 is applied through the entire RSG evolution phase.

Appendix D: Dust model interpolation

The interpolation of the DUSTY+MARCS (DM) models requires some extra details. First, all the MARCS models are downscaled to a spectral resolution of 1000. Then, all the DM models are calculated at the reference luminosity of L₀ = 1 L_⊙ and are later rescaled to the nominal luminosity of the RSG. Calculating all DM models at the same luminosity has the advantage that the interpolation between the models can be done directly. We interpolate the absolute magnitudes M_X in filters V and K, and the inner dust radius r₁, because all are later needed for calculating the RSG light curve. Each DM model (i.e., each combination of T_eff and log g) consists of 31 dust models with varying amounts of dust, represented by the optical depth τ linearly increasing from log τ = −1 to log τ = 2.

Fig. D.1. Kiel diagram showing the RSG pulsation cycle and the MARCS model coverage. Dots along the RSG track are spaced equally in time every 1/20 of the pulsation period. The crosses show the combinations of Teff and log g for which there are MARCS models available and for which we can compute DUSTY+MARCS (DM) models. After extending the grid by one grid spacing in each direction (green squares), we can cover 80% of the RSG pulsation cycle by linearly interpolating the newly obtained grid of DM models (shaded region).

Fig. D.2. Absolute V-band magnitude and the inner dust radius r1 for the DM models with log g = 0.0 as a function of effective temperature for various optical depths. For comparison, we extrapolate the DM models from log g = 1.0 and 2.0 to log g = 0.0 wherever the models are available.

Next, we extrapolate the DM models by one grid spacing in T_eff at constant log g, both at the low and at the high-temperature edge of the grid. Afterward, we extrapolate the DM models by one grid spacing in log g at constant T_eff toward the low surface gravity edge of the model grid. Fig. D.2 shows that the extrapolation along log g works well. Once we extend the grid, we can cover 80% of the RSG pulsation cycle with the DM model grid (Fig. D.1). The interpolation of the DM models is done linearly on the extended DM model grid.

To adjust the interpolated DM models to the nominal luminosity L of the RSG, we need to rescale the inner dust radius r₁ and absolute magnitudes M_X via

$$\begin{array}{cc}\hfill M_X& \to M_X-2.5{log}_{10}(L/L_0),\hfill \end{array}$$(D.1)

$$\begin{array}{cc}\hfill r_1& \to r_1\sqrt{L/L_0},\hfill \end{array}$$(D.2)

with the reference luminosity L₀ = 1 L_⊙ (Ivezić & Elitzur 1997).

To construct the light curve, we need to know the time evolution of the optical depth τ. The optical depth is not expected to be constant because the pulsating RSG changes both in luminosity and effective temperature. Hence, the inner dust radius r₁, defined at a dust condensation temperature of 800 K, needs to vary with time as well. First, we fix the optical depth τ₀ at maximum extent (ϕ = 0). For each DM model, we linearly interpolate along the τ axis to obtain the quantities at the desired optical depth. Based on τ₀, we can determine the dust density, ρ_x, at an arbitrary but fixed distance, x, which satisfies r₁ ≤ x ≤ 10³r₁ for all DM models simultaneously, via

$$\begin{array}{c}\hfill {\rho }_x={\rho }_1{\left(\frac{r_1}{x}\right)}^2.\end{array}$$(D.3)

The density ρ₁ at distance r₁ is given by

$$\begin{array}{c}\hfill {\rho }_1=\frac{\tau {\eta }_1}{r_1\mu /\rho }\propto \frac{\tau }{r_1},\end{array}$$(D.4)

with η₁ the normalized dimensionless density at r₁ (see Ivezić & Elitzur 1997, equation 4 for the definition of η) and μ/ρ the mass attenuation coefficient. Both η₁ and μ/ρ are constant amongst all DM models because we assume the same density scaling and dust composition. The power-law scaling of the density profile is set by the input physics of the DUSTY models (Sect. 2.3).

To obtain the τ evolution of the light curve, we require the density ρ_x to be constant throughout the pulsation cycle. With this assumption and together with Eq. (D.3), we can determine the optical depth τ for any DM model at given phase ϕ. The absolute magnitudes M_X are obtained by linear interpolation of the DM model along the τ axis. This way, M_X can be determined for all the points during the pulsation cycle, allowing us to construct the entire light curve.

Appendix E: Envelope structure during the pulsation cycle

We show various structure evolution diagrams of the envelope throughout the pulsation cycle in Figs. E.1 and E.2.

From the released recombination energy (Fig E.1a), one can see the different phases of recombination and re-ionization of hydrogen and helium (see also Figs E.1b-d). We find a short phase of recombination of He²⁺ at ϕ = 0.3 that leads to the expansion of intermediate layers. However, this expansion phase is quenched as the entire envelope above falls back. It seems that the natural oscillation frequency of an envelope driven by He²⁺ recombination is shorter than the pulsation period we find for the entire RSG envelope.

Fig. E.1. Structure evolution diagrams of the envelope of 2 full pulsation cycles. Colors show the specific released recombination energy ϵrec in panel (a), the ionization fractions of hydrogen and helium in panels (b)-(d), the radiative luminosity Lrad in units of the Eddington luminosity LEdd in panel (e), the opacity κ in panel (f), the specific entropy s in units of Avogadro’s number NA and the Boltzmann constant kB in panel (g), and the radial velocity v in units of the sound speed cs in panel (h). The dashed lines represent lines of constant mass for masses of 12, 11, 10, 9, and 8 M⊙ from top to bottom. For better visualization, the logmod transformation is used for the specific recombination energy, with logmod(x) = sign(x)log10(|x|+1) (John & Draper 1980).

Fig. E.2. Structure evolution diagrams similar to Fig. E.1. Colors show the temperature T in panel (a), the first adiabatic index Γ1 = (∂lnP/∂lnρ)s in panel (b), the density ρ in panel (c), the absolute radial velocity in units of the convective velocity vconv in panel (d), the total pressure P in panel (e), the radial velocity in panel (f), the fraction of the radiation pressure Prad of the total pressure in panel (g), and the velocity divergence ∇ ⋅ v in panel (h).

The neutral envelope falls back at highly supersonic speeds (Fig. E.1h) and is re-ionized at very large rates (≈10⁸ erg s⁻¹ g⁻¹). This causes the formation of shock waves (see Fig. E.2h that traces the shocks via the velocity divergence). In these layers, the opacity reaches the highest values throughout the envelope (Fig. E.1f). This also causes the radiative luminosity to locally exceed the Eddington luminosity (Fig. E.1). The re-ionization layer causes a sharp increase in density (Fig. E.2c) and pressure (Fig. E.2e). Ultimately, all these factors lead to the bounce-off and re-expansion of the outer layers, assisted by recombination. The constant-mass lines in Fig. E.1 clearly show that the layers further inside the envelope keep contracting. As they become more compressed, the temperature (Fig. E.2a) and the pressure increase, slowing down the layer’s contraction and starting its re-expansion at ϕ = 0.5. This turnaround process deep inside the envelope is much more gradual compared to the outer layers.

Around ϕ = 0.7, the outer layers of the envelope expand at supersonic velocities (Fig. E.1h). In some of these layers, energy transport is moderated via convection. We find that the expansion velocity of these layers exceeds the convective velocity that is predicted by mixing length theory (Fig E.2d). The application of mixing length theory to these layers is most certainly not appropriate, as we deviate far from hydrostatic equilibrium, and the simulation timescale is comparable to the convective turnover timescale.

Appendix F: Photosphere evolution

We show the time evolution of the photosphere mass coordinate m_phot, the photosphere radius r_phot, and the effective temperature T_phot in Fig. F.1. For comparison, we also show the bolometric light curves and the photospheric velocities (see Fig. 7).

For 20 d < t < 80 d, we see a clear ordering of m_phot for any fixed time. For a larger progenitor radius, we find a larger photosphere mass coordinate. This means that the photosphere recedes later for the more extended progenitors (points A, B, and F) compared to the compact progenitors (points C, D, and E). The rate at which the photosphere moves to lower masses is similar.

Fig. F.1. Time evolution of the photosphere mass coordinate (panel a), the radius of the photosphere (panel b), the effective temperature of the photosphere (panel c), and the photospheric velocity (panel d). The light dashed lines indicate the phase when the photosphere temperature drops below 103.75 K. At this point, SNEC has difficulties estimating the exact location of the photosphere (Morozova et al. 2015). The dash-dotted line illustrates the hydrostatic progenitor. The inset shows the radius evolution of one pulsation cycle and indicates the points A–F.

The photosphere radius of the most extended progenitor (point A) is the smallest until ∼30 d, but then expands to the largest overall radius. This is a direct consequence of the flat photospheric velocity evolution and a lower peak velocity compared to the more compact progenitors.

The SN from the most extended progenitors also shows the highest effective temperatures during the first 30 d. This is expected from the Stefan-Boltzmann law. The small photospheric radius and the large bolometric luminosity necessarily lead to a high effective temperature.

All Figures

Fig. 1. Radial pulsation of the RSG. The radius evolution is shown in panel a with a zoom-in on one pulsation cycle in panel b. The dash-dotted black line indicates the radius of the hydrostatic model. Panel c shows one pulsation cycle in the HRD. The arrow indicates the direction of the loop in the HRD, and the markers are spaced equally in time every 1/20 of the pulsation period. The hydrostatic model, the time-averaged effective temperature and luminosity, and their uncertainties are shown as well. For comparison, we show the typical uncertainty of the luminosity and effective temperature of Type II SN progenitors from Smartt (2015) as a gray marker. The colored markers labeled A–F in panels b and c indicate characteristic stages during the pulsation at which we calculate SN light curves. The pulsation phase, ϕ, is defined to be zero at maximum expansion.

In the text

Fig. 2. Restructuring of the envelope during the initial transient. Panel a shows the radius of the envelope (similar to Fig. 1a). The partial ionization zones and the mass coordinates of the recombination fronts of hydrogen and helium are shown in panel b, as well as the total mass of the star. The recombination fronts are defined where the ionization fraction of the respective species reaches 0.5. The partial ionization zones are defined to be the layers where the ionization fraction of the respective species is between 10 and 90%. Panel c shows the specific entropy in units of Avogadro’s number NA and the Boltzmann constant kB, as well as the density. Both quantities are reported at a mass coordinate of 5.5 M⊙, which is close to the bottom of the convective envelope. Panel d shows the radiative cooling timescale in units of the dynamical timescale. The total energy of the convective envelope is shown in panel e. The gray band indicates the catastrophic cooling event at t ≈ 32.5 yr.

In the text

	Fig. 3. Specific entropy, s, in units of Avogadro’s number, NA, and the Boltzmann constant, kB, at various mass coordinates inside the convective envelope. The gray band indicates the restructuring event at t ≈ 32.5 yr.
In the text

Fig. 4. Recombination structure of the RSG envelope at points A–F. The six colored circles show the radius of the RSG at points A–F. The coloring inside the circles shows the specific released recombination energy ϵrec as logmod(ϵrec/103 erg s−1 g−1) on a linear radial scale, with logmod(x) = sign(x)log10(|x|+1) (John & Draper 1980). Red corresponds to energy released by recombination, while blue corresponds to energy being removed via ionization. The square insets at each point show the recombination luminosity Lrec, i = ∫Menvϵrec, i dm for the species i = H+,  He+ and He2+ on a linear scale. Hatched regions indicate the partial ionization zone for hydrogen and helium, defined where the ionization fraction of the corresponding species is between 1% and 99%. Light gray shading at points A and B corresponds to neutral layers (ionization fractions less than 1%). The central plot shows the radius evolution of one pulsation cycle and indicates the 6 highlighted points A–F. The temporal evolution of the recombination structure is available as an online movie.

In the text

Fig. 5. Density profiles of the RSG envelope at points A–F during the pulsation cycle in terms of mass coordinate (panel a) and radius coordinate (panel b). The highlighted regions show hydrogen and helium partial ionization zones where the ionization fraction is between 10% and 90%. The density profile of the hydrostatic model is shown for comparison. Arrows along the profiles indicate the radial velocity, both the direction (inwards or negative; outwards or positive) and the magnitude, and highlight layers that are expanding/contracting. The inset in panel a shows the density structure of the interior layers, up to where the core was cut out. The shading indicates the core region of the RSG. The inset in panel b shows the radius evolution of one pulsation cycle and indicates the points A–F. The temporal evolutions are available as an online movie panel a and an online movie panel b.

In the text

Fig. 6. V- and K-band light curves of the pulsating RSG surrounded by a dust shell. The parameter τ0 sets the optical depth at the beginning of the pulsation cycle (ϕ = 0). The light dotted lines indicate regions outside the MARCS grid (too low Teff) that do not allow us to make photometric estimates. For better visualization, two full pulsation cycles are shown. The light, horizontal lines show the time-averaged magnitudes. For comparison, the phase-folded, visual light curve of VX Sgr is shown as well, with data taken between JD = 2640000 and 2644000, and assuming a period of 747 days.

In the text

Fig. 7. Bolometric light curves (panel a) and photospheric velocities (panel b) for SNe at points A–F during the pulsation cycle. The explosion of the hydrostatic model is shown for comparison. The inset shows the radius evolution of one pulsation cycle and indicates the points A–F. Once the photospheric temperature drops below 103.75 K, the velocities are no longer reliable because SNEC cannot accurately estimate the location of the photosphere due to limitations in low-temperature opacities (Morozova et al. 2015). The complete ϕ-evolution of the light curve and the photometric velocity is available as an online movie.

In the text

Fig. 8. Connecting feature in the SN light curves to features in the density profiles of the RSG envelope. Panel a is similar to Fig. 7a and panel b is similar to Fig. 5. The density profiles are shown as a function of the mass until the surface. Colored markers connect features in the density profiles to features in the SN light curves. Features in the light curve arise when the SN photosphere reaches the mass coordinates highlighted in the density profile.

In the text

	Fig. A.1. Comparison of the growth of the pulsations (panel a) and the steady state pulsations (panel b) between the undamped and damped model. The shading during the initial ∼10 yr indicates the magnitude of the damping factor γ.
In the text

Fig. A.2. Restructuring of the envelope for the model with initial damping, similar to Fig. 2. Panel (a) shows the specific entropy in units of Avogadro’s number NA and the Boltzmann constant kB, as well as the density. Both quantities are reported at a mass coordinate of 5.5 M⊙, which is close to the bottom of the convective envelope. The total energy of the convective envelope is shown in panel (b). For comparison, the undamped model is shown as well (same as Fig. 2c and e).

In the text

	Fig. B.1. Comparison of the growth of the pulsations (panel a) and the steady state pulsations (panel b) between the model with different cut-out core masses. The total helium core mass is 5.1 M⊙, defined by a hydrogen abundance less than 0.1.
In the text

	Fig. C.1. Comparison of the pulsations of two RSG models with different αMLT values over two pulsation cycles.
In the text

	Fig. C.2. Similar to Fig. 7a, but for the progenitor with αMLT = 3, where different colors correspond to different explosion phases. The light curves for the progenitor with αMLT = 1.8 are shown with gray colors as a comparison.
In the text

Fig. D.1. Kiel diagram showing the RSG pulsation cycle and the MARCS model coverage. Dots along the RSG track are spaced equally in time every 1/20 of the pulsation period. The crosses show the combinations of Teff and log g for which there are MARCS models available and for which we can compute DUSTY+MARCS (DM) models. After extending the grid by one grid spacing in each direction (green squares), we can cover 80% of the RSG pulsation cycle by linearly interpolating the newly obtained grid of DM models (shaded region).

In the text

	Fig. D.2. Absolute V-band magnitude and the inner dust radius r1 for the DM models with log g = 0.0 as a function of effective temperature for various optical depths. For comparison, we extrapolate the DM models from log g = 1.0 and 2.0 to log g = 0.0 wherever the models are available.
In the text

Fig. E.1. Structure evolution diagrams of the envelope of 2 full pulsation cycles. Colors show the specific released recombination energy ϵrec in panel (a), the ionization fractions of hydrogen and helium in panels (b)-(d), the radiative luminosity Lrad in units of the Eddington luminosity LEdd in panel (e), the opacity κ in panel (f), the specific entropy s in units of Avogadro’s number NA and the Boltzmann constant kB in panel (g), and the radial velocity v in units of the sound speed cs in panel (h). The dashed lines represent lines of constant mass for masses of 12, 11, 10, 9, and 8 M⊙ from top to bottom. For better visualization, the logmod transformation is used for the specific recombination energy, with logmod(x) = sign(x)log10(|x|+1) (John & Draper 1980).

In the text

Fig. E.2. Structure evolution diagrams similar to Fig. E.1. Colors show the temperature T in panel (a), the first adiabatic index Γ1 = (∂lnP/∂lnρ)s in panel (b), the density ρ in panel (c), the absolute radial velocity in units of the convective velocity vconv in panel (d), the total pressure P in panel (e), the radial velocity in panel (f), the fraction of the radiation pressure Prad of the total pressure in panel (g), and the velocity divergence ∇ ⋅ v in panel (h).

In the text

Fig. F.1. Time evolution of the photosphere mass coordinate (panel a), the radius of the photosphere (panel b), the effective temperature of the photosphere (panel c), and the photospheric velocity (panel d). The light dashed lines indicate the phase when the photosphere temperature drops below 103.75 K. At this point, SNEC has difficulties estimating the exact location of the photosphere (Morozova et al. 2015). The dash-dotted line illustrates the hydrostatic progenitor. The inset shows the radius evolution of one pulsation cycle and indicates the points A–F.

In the text