Instabilities in the Yellow Hypergiant domain

Yellow Hypergiants (YHGs) are massive stars that are commonly interpreted to be in a post-red supergiant evolutionary state. These objects can undergo outbursts on timescales of decades, which are suspected to be due to instabilities in the envelope. To test this conjecture, the stability of envelope models for YHGs with respect to infinitesimal, radial perturbations is investigated. Violent strange mode instabilities with growth rates in the dynamical regime are identified if the luminosity to mass ratio exceeds $\approx 10^4$ in solar units. For the observed parameters of YHGs we thus predict instability. The strange mode instabilities persist over the entire range of effective temperatures from red to blue supergiants. Due to short thermal timescales and dominant radiation pressure in the envelopes of YHGs, a nonadiabatic stability analysis is mandatory and an adiabatic analysis being the basis of the common perception is irrelevant. Contrary to the prevailing opinion, the models considered here do not exhibit any adiabatic instability.


INTRODUCTION
The upper domain of the Hertzsprung-Russell (HR) diagram is populated by massive stars (> 8 M ⊙ ) in their diverse evolutionary states.One of these stages comprise the yellow hypergiants (YHGs).These stars spread in temperature from about 4000 − 8000 K and have luminosities in the range 5.4 ≤ log L/L ⊙ ≤ 5.8.YHGs are commonly interpreted as being in their blue-ward evolution after having passed through a previous red-supgergiant (RSG) phase (e.g., de Jager 1998;Gordon & Humphreys 2019).According to stellar evolutionary tracks for rotating single stars (e.g. by Ekström et al. 2012), they might have evolved from progenitors with initial masses 20 − 40 M ⊙ because stars in this mass range are suggested to evolve back to the blue, hot side of the HR diagram where they may evolve into hot supergiants, such as the luminous blue variables, B[e] supergiants, or Wolf-Rayet stars.This prediction of the evolutionary models is consistent with the observational finding of an apparent absence of type II-P supernovae for stars that are initially more massive than ∼ 20 M⊙ and whose progenitors are RSGs (Smartt 2009).These evolutionary models also predict that stars more massive than 40 M ⊙ apparently do not evolve into cool RSGs.These stars seem to reach their turning point at significantly higher temperatures from where they evolve back to the blue side, and they do so the earlier (i.e. the hotter) the more massive they are.Consequently, evolution into a YHG is restricted to stars within a very narrow initial mass range, con-⋆ E-mail: wglatze@astro.physik.uni-goettingen.de(WG) sistent with the position of currently confirmed YHGs in the HR diagram (e.g., Kourniotis et al. 2022).
The distinction between pre-and post-RSG stars in the yellow domain is based on the significant spectroscopic and photometric variability of the latter.YHGs can also display indication for large-scale nebulae formed from the extensive massloss during the RSG phase.This is particularly true for the more massive objects, which transit much faster from the red to the yellow domain, leaving not much time for the released mass to expand and dilute.And indeed, the more massive stars in the currently known sample of Galactic YHGs display either an extended nebulocity (e.g., IRC +10420, Tiffany et al. 2010) or indication for large spherical shells of expanding cold molecular gas (e.g., IRC +10420 and HD 179821 Oudmaijer et al. 2009).
Besides these large-scale ejecta, YHGs can also be embedded in material that the stars have released more recently, most likely during one or more outburst events.Clear indication for such recent ejection episodes has been found for only a few cases, such as the Galactic object IRAS 17163-3907 (also known as the Fried Egg nebula, Lagadec et al. 2011) which seems to have experienced at least three outbursts within the past 100 years.These outbursts have led to the formation of three individual dust shells around the star (Koumpia et al. 2020).
But not every outburst releases enough mass to guarantee the production of significant amounts of detectable dust.The prime example is the Galactic object ρ Cas that has experienced at least four outbursts during the past ∼ 90 yr (Maravelias & Kraus 2022), but only after its outburst in 1946-47 emission from dust could be detected that must have formed from the released matter (Jura & Kleinmann 1990).Since then, this dust shell expanded and cooled, but no new dust has formed in detectable amounts from the more recent events (Shenoy et al. 2016).
Ejected circumstellar material can also be traced by static nebular line emission of low-excitation metal lines, such as Fe i, Sr ii, Y ii, and Ba ii (e.g., Lobel et al. 1998;Kourniotis et al. 2022), by emission of forbidden lines such as [Ca ii] and [O i], whereby [O i] is typically seen in hotter YHGs such as IRC +10420 and V509 Cas (Aret et al. 2017;Klochkova 2019), or by emission of warm molecular gas.The CO rovibrational bands are the most prominent molecular emission features.These molecular bands have been detected in the near-infrared spectra of a number of objects such as the Galactic YHGs V509 Cas and ρ Cas (Lambert et al. 1981;Gorlova et al. 2006;Kraus et al. 2022a), HD 179821 (Hrivnak et al. 1994), [FMR2006] F15 (Davies et al. 2008;Kraus et al. 2023) and two objects in the Large Magellanic Cloud (LMC), the stars HD 269723 and HD 269953 (McGregor et al. 1988;Oksala et al. 2013).Moreover, hot water vapor emission has recently been discovered from the environment of HD 269953 (Kraus et al. 2022b).
Outbursts in YHGs are usually identified by a sudden and steep decrease in visual brightness of the star along with indications for a rapid drop in spectroscopic temperature and the formation of TiO molecules in the expanding envelopes, which give rise to characteristic absorption bands in the spectra.The decrease in spectroscopic temperature makes the star to seemingly undergo a red-ward excursion in the HR diagram, until the episode of strong mass loss cedes and the material expands and dilutes.The recovery phase is usually much longer than the onset of the outburst and proceeds with a gradual brightening, back to the object's pre-outburst magnitude, and an apparent heating-up of the star causing its movement back to the hotter pre-outburst position in the HR diagram.The outburst duration of individual YHGs can be very different.Recorded have been events that lasted for decades as for example experienced by Var A in M33 (Humphreys et al. 2006) or for just a couple of years as for the Galactic star ρ Cas (Lobel et al. 2003;Kraus et al. 2019;Maravelias & Kraus 2022).
The cause of the outburst activity of YHGs has been ascribed to in the literature as due to instabilities occurring in the envelope or atmosphere of the stars (e.g., Nieuwenhuijzen & de Jager 1995;de Jager et al. 2001).In particular, it has been proposed that, when a YHG approaches a temperature of about ∼ 7000 K, its atmosphere starts to become unstable leading to substantial mass loss (e.g., Stothers & Chin 1993, 2001;de Jager et al. 2001;Lobel 2001).This temperature has been suggested to mark the lower boundary of a domain that has been dubbed the "yellow evolutionary void" (Nieuwenhuijzen & de Jager 1995;de Jager & Nieuwenhuijzen 1997) because of the apparent lack of stars observable within this region.Furthermore, the outburst activity of YHGs has been referred to as bouncing at the yellow void (de Jager & Nieuwenhuijzen 1997) or respectively the yellow/white wall (Oudmaijer & de Wit 2013), because of the apparent redward directed excursion the star undertakes after each event.
In the current work, we critically review the concept of the dynamical instabilities, in particular the adiabatic instability that is usually claimed to be responsible for the outbursts and the mass ejections in YHGs.On the basis of an adia-batic stability analysis we prove that all stellar models in the YHG domain are stable, questioning the existence of the yellow evolutionary void.Instead, we propose that the outbursts of YHGs could be related to strange-mode instabilities.As has been shown by Gautschy & Glatzel (1990b) and Glatzel (1994), the excitation of these modes is to be expected in massive stars with high values of their luminosity over mass ratio, for which post-RSGs are excellent candidates (Saio et al. 2013).Because strange-mode instabilities have the potential to trigger time-variable mass loss and mass ejections (Glatzel et al. 1999), they provide an excellent, alternative mechanism to drive outbursts in YHGs.We present a thorough stability analysis with respect to linear nonadiabatic radial perturbations focusing on the parameter space occupied by confirmed YHGs.We demonstrate the occurrence of strange-modes in all suitable models and show that their appearance does not (or only very mildly) depend on the effective temperature of the star.
In Section 2 the construction of stellar models and the basic assumptions and methods of the stability analysis are described.The results for models of ρ Cas and Yellow Hypergiants in the LMC together with an investigation of their dependence on effective temperature are presented in Section 3. Section 4 contains an extensive critical discussion in particular with respect to the common perception of adiabatic instability in the upper domain of the HR diagram.Our conclusions follow in Section 5.

Stellar models
In order to represent the observed properties of YHGs as accurately as possible, the study is based on envelope models constructed for the observed stellar parameters luminosity, effective temperature and chemical composition.The uncertainty in mass is taken into account by considering wide mass ranges which should include the values suggested by both the spectroscopic analysis and the comparison with evolutionary models.In order to demonstrate the dependence on effective temperature of the results of the stability analyses, we shall also consider model sequences with varying effective temperature and fixed luminosity, chemical composition and mass.
For prescribed luminosity, effective temperature, chemical composition and mass the structure of the stellar envelope between the photosphere and some suitably chosen bottom boundary can be determined by initial value integration of the equations of hydrostatic equilibrium, energy transport and mass conservation, where unambiguous initial values are imposed at the photosphere.By definition the luminosity is constant throughout the envelope.The bottom boundary is defined in terms of a maximum temperature which guarantees that nuclear burning does not prevail.It corresponds to a finite radius.
Concerning the treatment of convection, its onset is determined by Schwarzschild´s criterion, standard mixing length theory (Böhm-Vitense 1958) with 1.5 pressure scale heights for the mixing length is adopted for its description, and overshooting as well as semiconvection are disregarded.Opacities have been taken from the OPAL tables (see Rogers & Iglesias 1992, Iglesias & Rogers 1996and Rogers et al. 1996).

Stability analysis
In the present study, we test the envelope models for YHGs for stability with respect to infinitesimal radial perturbations.The associated mathematical problem is derived and described, e.g., in Baker & Kippenhahn (1962).Adopting their notation and treating convection according to the "frozen -in -approximation" (see, e.g., Baker & Kippenhahn 1965), the boundary eigenvalue problem posed by the analysis of radial linear nonadiabatic stellar stability and pulsations (LNA analysis) is solved using the Riccati method (see Gautschy & Glatzel 1990a).In addition to the LNA analysis, the envelope models have been subject to a standard radial linear adiabatic stability analysis (for details see Cox 1980).
As a result of the stability analyses, we obtain for each envelope model its complex eigenfrequencies σ, where the real parts σr correspond to the pulsation frequencies, and the imaginary parts σi indicate the growth or damping rates of the various modes.In our normalisation σi > 0 corresponds to damping (stability), σi < 0 to growth and excitation (instability).For convenience, the eigenfrequencies will be presented in dimensionless form, i.e., they will be normalised by the global free fall time (cf.Baker &Kippenhahn 1962 andGautschy &Glatzel 1990a).This normalisation is common for theoretical studies such as the present investigation.In particular, it avoids the masking of results by the period density relation.
In the radial linear adiabatic analysis (also referred to as the adiabatic approximation), the boundary eigenvalue problem is equivalent to a (selfadjoint) Sturm -Liouville problem (see, e.g., Ledoux & Walraven 1958or Cox 1980).As a consequence, σ 2 is real and forms an infinite, well ordered sequence with a smallest (fundamental) eigenvalue and a limit point at infinity in the adiabatic approximation.Thus σ 2 < 0 for the fundamental eigenfrequency is a necessary and sufficient condition for instability (and vice versa) in the adiabatic approximation.Accordingly, adiabatic stability and instability can be determined by merely considering the fundamental adiabatic eigenfrequency.Instability sets in through σ = 0.

ρ Cas
Adopting observed values for the luminosity (L = 5 • 10 5 L⊙, Humphreys 1978), the mean spectroscopic effective temperature (T eff = 7000K, Lobel et al. 1994;Kraus et al. 2019), solar chemical composition ((X, Y, Z) = (0.74, 0.24, 0.02)), and a range in mass between 19M⊙ and 50M⊙, including the most likely value of the star's current evolutionary mass of 24.1M⊙ (Kraus et al. 2019), a sequence of envelope models with the mass as a parameter has been constructed and tested for stability.Real and imaginary parts of the lowest order eigenfrequencies σ are presented as a function of mass in Figs. 1 and 2.
At high masses, all modes are damped and their frequencies are regularly spaced, as expected for an ordinary acoustic resonator.With decreasing mass (below ≈ 35M⊙), multiple mode crossings and mode pairings unfolded both according to the ´avoided crossing´and the ´instability bandć oupling schemes (cf., e.g., Gautschy & Glatzel 1990b) are found, which are associated with the occurrence of instabili- ties having growth rates in the dynamical regime.Apart from one strongly damped mode, whose frequency and damping increases, frequencies and dampings tend to decrease with decreasing mass.For masses below ≈ 25M⊙, damped and unstable modes exhibit an approximately complex conjugate symmetry, which is typical for the pure form of mode coupling according to the ´instability band´scheme.
The behaviour of modes and the occurrence of instabilities is a consequence of the change with mass of the stellar structure.The latter is shown in Fig. 3 by means of the density stratification of envelope models for ρ Cas with four different masses.
The core -envelope structure of these models, where a small core contains almost the total mass of the star and a tenuous envelope with negligible contribution to the stellar mass covers the entire space, is becoming more and more pronounced as the mass decreases.This change in structure, in particular the decrease of density in the envelope with decreasing mass, has a direct impact on the contribution of gas pressure to total pressure in the stellar envelope.The ratio β of gas pressure to total pressure for the models considered is presented in Fig. 4.
Whereas for M = 50M⊙ the fraction of gas pressure is still higher than 25% everywhere, it falls below 10% in almost the entire envelope for M = 19.1M⊙.Thus the envelopes studied become more and more dominated by radiation pressure as the mass decreases.Another crucial consequence of the structure and low densities of the envelopes considered refers to the timescales governing acoustic waves in the stellar envelope (see also the analogue discussions in Gautschy &Glatzel 1990b andGlatzel 2021).
The local dynamical timescale may be estimated as the time needed by a sound wave to cross a mass shell with thickness ∆r.Estimating the sound speed as c 2 Sound ∝ p/ρ (p: pressure, ρ: density), it is given by: τDyn ∝ ∆r ρ/p . (1) On the other hand, the local thermal timescale of a mass shell with mass ∆m may be defined as the time needed to radiate its thermal energy content at the local luminosity, where the thermal energy content might be expressed as the product of the specific heat cp, the temperature T and the mass ∆m.Rewriting the latter in terms of the density ρ and the volume of the mass shell, we finally obtain for the local thermal time scale: (2) Both the local dynamical and the local thermal timescale depend on the thickness ∆r of the mass shell considered.Unless there are further arguments how to choose ∆r, they can be given any value since the choice of ∆r is ambiguous.Thus the local dynamical and thermal timescales given by equations ( 1) and ( 2) are ill defined quantities without any physical relevance.However, their ratio being independent of ∆r is well defined and given by: The ratio of thermal and dynamical timescales as a function of relative radius for the envelope models discussed is shown in Fig. 5.As in any stellar model, this ratio achieves very high values in the core and is smallest (maybe even below unity) at the surface.As a consequence, all sound waves become adiabatic in the deep interior of the star and nonadiabatic effects have to be taken into account in a certain domain below the stellar surface, where the ratio of thermal and dynamical timescales is of order unity or smaller.This domain depends on the stellar model and shrinks with increasing mass for the models discussed (see Fig. 5).In other words, with decreasing mass we expect the adiabatic approximation to become less and less valid.Instead of characterising them by an infinite thermal timescale (adiabatic approxiation), low mass models for ρ Cas should rather be described by the opposite approximation of a small or vanishing thermal timescale.The latter corresponds to the non-adiabatic reversible (NAR) approximation (see, e.g., Gautschy & Glatzel 1990b).
In the NAR approximation, eigenfrequencies occur in complex conjugate pairs, i.e., modes are either neutrally stable, or damped and unstable modes appear simultaneously thus forming pairs with the same frequency and identical moduli of growth and damping rates.Motivated by the considera- tion that the NAR approximation might be applicable to low mass models for ρ Cas, we have performed a corresponding analysis whose results are shown in Figs. 6 and 7.
In the NAR approximation, for high masses the modes are found to be neutrally stable.With decreasing mass below ≈ 25M⊙ some of the adjacent modes merge to form complex conjugate pairs of damped and unstable modes in a way that is quite close to the analysis without approximation (cf.Figs. 1 and 2).For comparison, Fig. 8 contains the imaginary parts of σ both from the exact and the NAR analysis.With respect to the strong dynamical instabilites at low masses,  and 7, but with results from the exact analysis (thin lines) and using the NAR approximation (thick lines) superimposed.Note the similarity of exact and NAR results, in particular in the domain of strong instabilites with growth rates in the dynamical regime.
their complex conjugate symmetry and the associated mode interactions, we conclude that the NAR approximation provides at least qualitatively satisfactory results.These findings thus also support the assumption that the opposite approximation of an infinite thermal timescale, the adiabatic approximation, should be invalid for the models considered.
To prove this conjecture, an adiabatic analysis has been performed for the ρ Cas models.The results in terms of the frequencies of the three lowest neutrally stable adiabatic modes are shown and compared with the exact frequencies in Fig. 9. From Fig. 9 we deduce that -as expected -the adiabatic frequencies do not provide an approximation to the exact frequencies in any respect, not even for high masses. .Same as Fig. 1, but with the three lowest neutrally stable adiabatic eigenfrequencies (dashed lines) added.We emphasize that there is no instability in the adiabatic approximation and the adiabatic frequencies do not appear to provide an approximation to the correct frequencies at all.
Moreover, the fundamental adiabatic mode which indicates adiabatic stability and instability, respectively, exhibits finite frequency for all models and thus no evidence at all for adiabatic instability.
To complete the discussion of the adiabatic analysis, the adiabatic exponent γ ad is shown for four ρ Cas models with different mass in Fig. 10.All models exhibit zones in which γ ad falls below the critical value 4/3.They are associated with the various ionization processes, each of them causing a minimum of γ ad .Even if these zones with γ ad < 4/3 do exist, their strength is not sufficient for adiabatic instability, i.e., the pressure weighted volumetric mean of γ ad does not fall below 4/3 (which would be sufficient for adiabatic instability).The fact that the fundamental mode has not shown any evidence for instability explicitly proves that the correct mean of γ ad is bigger than 4/3.Outside the ionization zones, γ ad is bigger than but close to 4/3 and increases with mass.This is a consequence of dominant radiation pressure (cf.Fig. 4): The limit of pure radiation pressure (β → 0) implies γ ad → 4/3.With increasing mass, β increases (see Fig. 4) and, together with it, also γ ad .

Yellow Hypergiants in the LMC
In this section, the dependence on metallicity of instabilities in the YHG domain will be studied.This is motivated by the recent investigations of Kourniotis et al. (2022) of evolved hypergiants in the LMC who classified the star HD 269723 as a luminous post-RSG and HD 271182 as a ρ Cas analogue and thus both as YHGs.These objects hence serve as ideal targets for our analysis.Both stars have a similar temperature of T eff ∼ 6000K but different luminosities.Kourniotis et al. (2022)  Figs.11 and 12 may be compared with their counterparts for ρ Cas, Figs. 1 and 2. Qualitatively, there is no difference between the results for ρ Cas and the LMC object.Mode interactions and associated instabilities do occur in the same way for both sequences.Instability sets in below ≈ 30M⊙ for the LMC object, at a somewhat lower mass than for ρ Cas, which is likely to be due to its smaller luminosity.As for ρ Cas, we have performed an adiabatic analysis for the LMC models.Its result for the sequence with L = 4.5 • 10 5 L⊙ in terms of the three lowest neutrally stable adiabatic eigenfrequencies is shown and compared with the exact results in Fig. 13 (cf.the counterpart for ρ Cas, Fig. 9).Except for the fundamental adiabatic mode at high masses, the adiabatic frequencies do not provide an approximation to the correct frequencies.Moreover, we emphasize that an adiabatic instability does not exist.
The results obtained for the sequence with L = 6•10 5 L⊙ are very similar to those for its counterpart with L = 4.5 • 10 5 L⊙, such that a separate discussion is redundant.Accordingly, for this sequence we only show the imaginary parts of σ as a function of mass in Fig. 14.As a consequence of the higher luminosity, the upper limit in mass for instability has shifted to higher masses (to around ≈ 40M⊙, compare Figs. 14 and  12).

Dependence on effective temperature of instability
In this section, the dependence on effective temperature of the instabilities found in the YHG domain will be studied.For this purpose, sequences of envelope models with the effective temperature as a parameter covering the range between red and blue supergiants are constructed and tested for stability.Adopting the chemical composition (X, Y, Z) = (0.74, 0.24, 0.02), the various sequences are characterized by the values selected for mass and luminosity.For a sequence .Same as Fig. 11, but with the three lowest neutrally stable adiabatic eigenfrequencies (dashed lines) added.We emphasize that there is no instability in the adiabatic approximation and the adiabatic frequencies do not appear to provide an approximation to the correct frequencies at all. with luminosity L = 5 • 10 5 L⊙ and mass M = 25M⊙, real and imaginary parts of the lowest order eigenfrequencies σ are presented as a function of effective temperature in Figs.
15 and 16.In addition to the nonadiabatic LNA analysis, also an adiabatic analysis has been performed.Its result, i.e., the frequencies of the two lowest neutrally stable adiabatic eigenfrequencies is also shown in Fig. 15.Figs. 15 and 16 demonstrate that the mode interactions and associated instabilities with growth rates in the dynamical regime identified in the YHG domain persist for the entire effective temperature range from red to blue supergiants.Whether the stable gap around log T eff ≈ 3.7 is significant, remains to be seen.Again, the adiabatic eigenfrequencies do not provide an approximation to the exact frequencies in any respect, as the adiabatic approximation does not apply.We emphasize that for the entire range of effective temperatures studied adiabatic instability does not exist.The treatment of dominant convection, in particular the coupling of strong convection and pulsation is still an unsolved problem and becomes important at the low temperature end of the model sequence.With respect to these uncertainties, the results in the RSG domain should be interpreted with caution.
With regard to the controversial discussion of adiabatic instability of massive stars in the BSG and LBV phase (see Glatzel &Kiriakidis 1998 andStothers 1999), we have performed an adiabatic analysis for additional sequences of envelope models.For four sequences, the results in terms of the lowest order adiabatic eigenfrequencies as a function of effective temperature are shown in Figs. 17 and 18. σ 2 is real and remains positive in any case.Thus all modes are neutrally stable and an adiabatic instability does not exist.The sequences of avoided crossings appearing at high effective temperatures in Fig. 17 are caused by the crossing of two sets of acoustic modes associated with two acoustic cavities in the stellar envelope and have been discussed, e.g., by Kiriakidis et al. (1993) and Glatzel & Kiriakidis (1998).

Nonadiabatic stability analysis
The study of the stability of massive stars in the vicinity of the Humphreys -Davidson limit dates back to the investigation of Glatzel & Kiriakidis (1993), where violent instabilities with growth rates in the dynamical range have been identified.These instabilities have been found to be associated with mode coupling and the apparent occurrence of additional unexpected, until then incomprehensible modes, which were addressed as ´strange modes´, the associated instabilities as ´strange mode instabilities´.With the arrival of new opacities which correctly account for the contribution of heavy elements (Rogers & Iglesias 1992), the study by Glatzel & Kiriakidis (1993) has been repeated by Kiriakidis et al. (1993) on the basis of these opacities with similar results for effective temperatures below ≈ 20 000K.For higher effective temperatures additional metallicity dependent strange mode instabilities have been identified.The boundary in the HRD, above which all stellar models independent of metal- licity exhibit violent strange mode instabilities with growth rates in the dynamical regime, was found to coincide with the observed Humphreys -Davidson limit.Meanwhile the investigations have been confirmed several times and strange mode instabilites are in general expected to occur, if the luminosity to mass ratio (in solar units) exceeds values around ≈ 10 4 (see, e.g., Saio 2011).
According to the previous studies, strange mode instabilites are to be expected for the models of YHGs considered here, as their luminosity to mass ratio (in solar units) is of the order of ≈ 10 4 .In fact, the multiple mode couplings associated with violent instabilities which have been discussed in the previous sections correspond to strange modes and strange mode instabilities.That certain sequences of mode interactions may create the impression of additional ´strangeḿ odes, has been described by Gautschy & Glatzel (1990b) and may be seen in Fig. 15, where the eigenfrequencies are presented as a function of effective temperature rather than as function of mass.However, not all ´strange modes´are caused by mode coupling processes.E.g., the mode with the lowest frequency at M = 50M⊙ in Fig. 1 and strongly increasing frequency and damping rate for decreasing mass, belongs to the thermal rather than to the acoustic spectrum.This classification is suggested by considering the thermal timescale, being small for low masses and increasing with mass (cf.Fig. 5).Consequently, thermal frequencies should decrease with mass, as seen for the mode discussed.Note that around ≈ 35M⊙ the thermal mode interacts with the acoustic spectrum through avoided crossings, which for the lowest acoustic mode is not very well pronounced.
Another indication for the presence of thermal modes is found in Figs. 15 and 16.Around log T eff ≈ 4 both the real and imaginary part of the lowest order mode approach zero thus suggesting a classification as a thermal mode at least in a gap around log T eff ≈ 4. Across this gap the mode has not been followed continuously, since multiple interactions with thermal modes close to zero frequency make an unambiguous continuous tracking extremely difficult.Moreover, for the models considered the thermal spectrum is not decisive for stellar stability, even though its study may be of academic interest.

Adiabatic stability analysis
Due to short thermal timescales for large fractions of the considered stellar envelopes (see Fig. 5), the adiabatic analysis is invalid and has been found not to provide a satisfactory approximation to the exact results in any respect.Moreover, none of the models studied exhibits an adiabatic instability.As a value below 4/3 of the mean of the pressure -weighted volumetric mean of the adiabatic exponent is a sufficient condition for instability (see, Ledoux & Walraven 1958), this result strictly proves that the mean adiabatic index stays above 4/3 in any case.We note that the primary procedure to test adiabatic stability consists of calculating the (real) square σ 2 of the fundamental adiabatic eigenfrequency σ.Its sign provides a necessary and sufficient condition for adiabatic stability: For σ 2 < 0 we have instability, otherwise stability.
In the present study we have applied this method.A secondary procedure providing a sufficient condition for instability only is based on the Rayleigh -Ritz variational principle associated with the differential adiabatic perturbation problem (see Ledoux &Walraven 1958 andGlatzel &Kiriakidis 1998).It allows for an estimate of an upper bound for σ 2 of the adiabatic fundamental mode, whose sign is determined by 3⟨γ ad ⟩ − 4, where ⟨γ ad ⟩ denotes the pressure -weighted volumetric mean of γ ad over the entire stellar model.Thus, ⟨γ ad ⟩ < 4/3 is a sufficient but not necessarily necessary condition for instability.We emphasize that for the mean of the adiabatic exponent the entire stellar model has to be taken into account and that the domain for calculating the mean cannot be chosen arbitrarily.Otherwise ⟨γ ad ⟩ could be given any arbitrary value (see Fig. 10).
Our results concerning adiabatic stability are in blatant contradiction to the common perception of instability regions in the upper HR diagram (see, e.g., de Jager et al. 2001 andStothers &Chin 2001).The prevailing opinion assumes that two domains of adiabatic instability exist, where the ´blued omain accommodates the LBVs and the ´yellow -red´domain the YHGs.Apart from Stothers & Chin (2001) whose study reportedly is based on an explicit solution of the adiabatic wave equation, all other studies rely on a consideration of the mean adiabatic exponent for the determination of the domains of instability.
With respect to the ´blue´instability domain, we have already contradicted the common view in an earlier paper (Glatzel & Kiriakidis 1998) (a) by arguing that the adiabatic approximation is not valid at all due to short thermal timescales in the stellar envelope, and by (b) explicitly demonstrating that there is no adiabatic instability.Thus the arguments raised here are not new and have been produced already in earlier studies in a similar context.

The relation between nonadiabatic and adiabatic stability analysis
Concerning the objections raised in Glatzel & Kiriakidis (1998), comments were published by Stothers (1999), which also form the basis of the study by de Jager et al. (2001).
According to Stothers (1999), adiabatic stability and nonadiabatic stability are considered separately and independently.We strongly disagree, since the adiabatic analysis is an approximation and subordinate to the nonadiabatic analysis, which can only be applied if the thermal timescales tend to infinity.It depends on the stellar model, whether a complete nonadiabatic analysis is required or the adiabatic analysis is sufficient.For massive stars having short thermal timescales in their envelopes, the adiabatic approximation is not valid and the appropriate analysis must take thermal effects into account.There is no choice concerning the analysis, and the two approaches (nonadiabatic and adiabatic analysis) cannot be applied independently and separately.In particular, it is meaningless to distinguish between dynamical and pulsational instability.Note that in our studies the term ´dynamical´only refers to the timescale and not to the type or mechanism of an instability.In the notation of Stothers (1999) it seems that the term ´dynamical instability´is also used as an equivalent for ´adiabatic instability´.Thus we conclude that an adiabatic stability analysis for the massive stars considered is irrelevant.Compared to this, the fact that we do not find adiabatic instabilities is of minor importance.An explanation for the discrepancy is not presented by Stothers (1999).

Criteria for adiabatic instability and their derivation
In general, the criterion for adiabatic instability involving the adiabatic exponent is derived from the Rayleigh -Ritz variational principle associated with the differential adiabatic perturbation problem (see Ledoux &Walraven 1958 andGlatzel &Kiriakidis 1998).Stothers (1999) presents an alternative derivation providing the same instability criterion as a result.It is based on an integral relation for the eigenfrequency (equation (2) of Stothers 1999), which may be derived either from the virial theorem or from an integration of the adiabatic wave equation (see Ledoux & Walraven 1958).Contrary to the Rayleigh -Ritz variational principle, this relation is not quadratic but linear in the Lagrangean displacement, which here is to be regarded as a solution of the wave equation rather than as a test function of the variational principle.A simple inspection of the wave equation shows, that for γ ad = 4/3 and σ 2 = 0 a constant Lagrangean displacement provides a solution of the perturbation problem.Stothers (1999) 3) are in fact identically satisfied.However, ignoring and relaxing the initial assumptions γ ad ≈ 4/3 and σ 2 ≈ 0 (which guarantee the Heaviside function as a solution of the wave equation), to derive a relation for γ ad ̸ = 4/3 and σ 2 ̸ = 0, implies a contradiction.Accordingly, even if the result appears to be correct, we consider its derivation to be wrong.

Comparison with observable quantities
For a comparison with observable quantities it might be tempting to convert the dimensionless frequencies derived here into pulsation periods and e-folding times.However, we emphasize that the present study is entirely based on a linear analysis which does not allow for a determination of amplitudes.Therefore, e-folding times can, in principle, not be related to any observed light variations, even if there might be a tendency that final amplitudes reached in the nonlinear regime increase with the growth rate of the underlying instability (see, e.g., Grott et al. 2005).Nonlinear simulations of strange mode instabilities have shown (see, e.g., Glatzel 2009) that due to their strength the stellar envelope is inflated by successive shock waves.As a consequence, the pulsation period increases and finally is significantly larger than the period determined by the linear analysis.We expect this behaviour also for the strange mode instabilities in YHGs, except possibly for models close to the onset of instability, where the growth rates are small.Accordingly, a comparison of linear periods with observed periods must be treated with caution.
Taking these cautionary remarks into account, the dimensionless frequencies σ can be converted into pulsation periods P and e-folding times τe by where R denotes the stellar radius and G is the gravitational constant.Equation (4) represents the period density relation and contains the global dynamical timescale which may be expressed as Using equations ( 4) and (6), we obtain from Fig. 1 theoretical pulsation periods for ρ Cas in the range between approximately 16 and 292 d.Noteworthy, the latter agrees considerably well with cyclic photometric variability of 200-300 d reported for ρ Cas over the past decades (Zsoldos & Percy 1991;Percy et al. 2000;Kraus et al. 2019).Moreover, the models for stars with LMC metallicity (Fig. 11) predict periods of about 750 d from the lowest order eigenfrequency for stars with about 30 M⊙, which is similar to the dominant period derived by Kourniotis et al. (2022) from the photometric light curves of the two LMC objects HD 269723 (800 d) and HD 271182 (833 d).

CONCLUSIONS
We have investigated the stability of stellar models in the Yellow Hypergiant domain with respect to infinitesimal radial perturbations.For luminosity to mass ratios above ≈ 10 4 violent strange mode instabilities with growth rates in the dynamical regime have been identified.Adopting the luminosities and masses derived from observations, we thus predict YHGs to suffer from these instabilities.For luminosity to mass ratios above ≈ 10 4 the strange mode instabilities persist over the entire range of effective temperatures from RSGs to BSGs, except possibly for a small stable gap around log T eff ≈ 3.7.Whether this gap is significant, remains to be studied.Should it be relevant for stellar evolution, it could mean that stars are forced to evolve into this gap, and may be pushed back into the gap once they try to evolve into the surrounding unstable domains.
The envelopes of YHGs with a pronounced core envelope structure are characterized by short thermal timescales and dominant radiation pressure, which according to a model for strange mode instabilites by Glatzel (1994) are essential ingredients for the occurrence of strange mode instabilities.In accordance with the short thermal timescales the NAR approximation (vanishing thermal timescale) has been found to describe mode interactions and instabilities at least qualitatively correct, when compared with the exact results.In contrast to the NAR approximation the opposite approximation of infinite thermal timescale (adiabatic approximation), as expected, does not provide an approximation to the exact results in any respect.We emphasize that due to the short thermal timescales in YHG envelopes a nonadiabatic analysis is mandatory and an adiabatic analysis is irrelevant.
According to the common perception, adiabatic instability causes instability regions in the upper HR diagram which also cover the YHG domain.Therefore, we have performed an adiabatic analysis for YHG models, even if the short thermal timescales indicate that the adiabatic approximation does not hold there.Contrary to the prevailing opinion, our results do not exhibit any adiabatic instability.Thus we disagree with the common conception in two respects: (a) The adiabatic approximation is not applicable, and (b) even if the adiabatic approximation was applicable, there is no adiabatic instability.
The linear analyses performed here do neither provide information on the amplitude that an unstable perturbation may reach, nor on the final fate of an unstable object.To determine them, the evolution of instabilities into the nonlinear regime needs to be followed by numerical simulation.For strange mode instabilites in massive stars such simulations (see, e.g., Glatzel et al. 1999 andGrott et al. 2005) indicate that pulsationally driven mass loss may be a consequence of the instability.Currently we are performing corresponding numerical simulations for YHG models, which we expect to provide information about whether strange mode instabilities can contribute to the observed mass loss and outbursts of these stars.Their results will be commented on in a subsequent paper.
Strange mode instabilities are not restricted to radial perturbations and occur also for nonradial perturbations in a similar way for low harmonic degrees l up to l ≈ 300 (see, e.g., Glatzel & Gautschy 1992, Glatzel & Mehren 1996and Glatzel & Kaltschmidt 2002).Thus we expect nonradial strange mode instabilities to be present also in YHG models.A corresponding study will be presented elsewhere.

Figure 1 .Figure 2 .
Figure 1.Real parts of the lowest order eigenfrequencies σ as a function of mass for envelope models with parameters resembling those of ρ Cas.Unstable modes are indicated by thick lines.

Figure 7 .Figure 8 .
Figure7.Same as Fig.2, but using the NAR approximation.Modes are either neutrally stable or appear in complex conjugate pairs of damped and unstable modes.
Figure9.Same as Fig.1, but with the three lowest neutrally stable adiabatic eigenfrequencies (dashed lines) added.We emphasize that there is no instability in the adiabatic approximation and the adiabatic frequencies do not appear to provide an approximation to the correct frequencies at all.
Figure13.Same as Fig.11, but with the three lowest neutrally stable adiabatic eigenfrequencies (dashed lines) added.We emphasize that there is no instability in the adiabatic approximation and the adiabatic frequencies do not appear to provide an approximation to the correct frequencies at all.

Figure 15 .Figure 16 .
Figure15.Real parts of the lowest order eigenfrequencies σ as a function of effective temperature for envelope models with luminosity L = 5 • 10 5 L ⊙ , mass M = 25M ⊙ and chemical composition (X, Y, Z) = (0.74, 0.24, 0.02).Unstable modes are indicated by thick lines.Dashed lines correspond to the two lowest neutrally stable adiabatic eigenfrequencies.

Figure 17 .
Figure17.The real and positive square σ 2 of the three lowest neutrally stable adiabatic eigenfrequencies as a function of effective temperature for envelope models with chemical composition (X, Y, Z) = (0.74, 0.24, 0.02), luminosity L = 10 6 L ⊙ and mass M = 70M ⊙ (full lines).Dashed lines represent the eigenfrequencies for envelope models with luminosity L = 5 • 10 5 L ⊙ and mass M = 45M ⊙ .

Figure 18 .
Figure18.The real and positive square σ 2 of the two lowest neutrally stable adiabatic eigenfrequencies as a function of effective temperature for envelope models with chemical composition (X, Y, Z) = (0.74, 0.24, 0.02), luminosity L = 10 6 L ⊙ and mass M = 40M ⊙ (full lines).Dashed lines represent the eigenfrequencies for envelope models with luminosity L = 5 • 10 5 L ⊙ and mass M = 25M ⊙ .